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Djelo Sage računalno okruženje za fizičare, čiji je autor Krešimir Kumerički, ustupljeno je pod međuna¬ 
rodnom licencom Creative Commons "Imenovanje-Dijeli pod istim uvjetima" (Attribution-ShareAlike) 
4.0. Za uvid u tu licencu, posjetite http : //creativecommons . org/licenses/by-sa/4.0/ 
deed. hr. 


Premda je ovaj dokument nastao radi potreba studija fizike na Fizičkom odsjeku Prirodoslovno- 
matematičkog fakulteta Sveučilišta u Zagrebu, največi dio (zapravo sve osim 5. poglavlja) bi trebao 
biti razumljiv i upotrebljiv svakome tko želi naučiti raditi u Sageu. 

Ispisi računalnog koda u ovom dokumentu su dobiveni unutar Sage računalnog okruženja, verzija 6.9. 
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Sadržaj 



POGLAVLJE 1 


Uvod 


1.1 Python vs. IPython vs. Sage 

• Python je računalni jezik opće namjene. Zahvaljujući, s jedne strane, dobroj čitljivosti koda, a 
s druge brojnim kvalitetnim bibliotekama za numeriku, simboliku, crtanje itd., Python se često 
koristi u istraživanjima i edukaciji iz fizike i srodnih područja znanosti i tehnike. 

• IPython je sučelje za Python namjenjeno interaktivnom radu. Postoje inačice za terminalski rad, 
grafičko sučelje i sučelje putem WWW preglednika. 

• Sage je matematički softver koji udružuje niz postojećih matematičkih (i drugih) biblioteka u 
zajedničko sučelje zasnovano na Pythonu, s ciljem kreiranja alternative komercijalnim softverima 
poput Mathematice, Maplea ili Matlaba. 

Zbog njegove univerzalnosti, u ovom dokumentu koristit ćemo uglavnom Sage. 


1.2 Zašto Python/Sage? 

1. Python je izrazito elegantan za upotrebu. Popularan je u znanstvenoj zajednici, pa postoji velik 
broj korisnih Python biblioteka i za specijalizirane znanstvene namjene.. 

2. Znanost mora biti reproducibilna i u načelu vječna. Računalni kod koji ovisi o softveru “zatvore¬ 
nog” koda to onemogućuje. 

3. Sage je najrazvijeniji sustav za računalnu algebru (CAS - Compute Algebra System) otvorenog 
koda {open source) 

Za detaljniju diskusiju o izboru softvera za računanje u znanosti, vidi Johansson, Introduction to scien- 
tific computing with Python, te, također, Koepke, 10 Reasons Python Rocks for Research (And a Few 
Reasons it Doesn’t). 
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POGLAVLJE 2 


Sučelje 


2.1 Sage radni list i ćelije 

Radni list {worksheet) je skup računalnog Sage koda s rezultatima i tekstualnim opisom koji je 

• prikazan na jednoj stranici u WWW pregledniku 

• može se pohraniti u jednu datoteku (ekstenzija . sws) 

Svi radni listovi jednog korisnika dostupni su putem poveznice Home na vrhu svakog radnog lista. 
Čitavo ovo grafičko sučelje naziva se Sage bilježnica {Sage notebook). (Postoji i obično tekstualno 
terminalsko sučelje.) 

Račun, tj. računalni kod radnog lista je organiziran u tzv. ćelije {cell). Svaka računska (input) ćelija je 
ograđena pravokutnikom i izvršava se kao cjelina na jedan od dva načina: 

1. pritiskom na kombinaciju tipki shift-Enter dok je kursor bilo gdje u ćeliji (Pritisak samo na 
Enter otvara novi redak u ćeliji), 

2. klikom mišem na gumb Evaluate koji se prikazuje ispod aktivne ćelije. 

Sage ispisuje rezultat neposredno ispod ćelije. Ako se eksplicitno ne zatraži drugačije (npr. naredbom 
print), bit će ispisan samo rezultat zadnjeg računa/komande/retka u toj ćeliji (ali izvršit će se naravno 
cijela ćelija). 




^ ^ 1 (ŠS . calculon.phy.hr:8080/home/kkumer/138/ 

*1 > ® H 


sDge. Tlie Sage Notebook kkumer Toegle | Home | Rublisbgd | Lag | Settmes | HdE I Rep.ort.a.Pr.Dbleni | Sigivo.ut 


Naslov radliog lista | Save |fŠave & guit I Tpiscard & guit | 

iFile... ^lAction... ’^IlData... I sage ^ - Typeset 


5-1 
1 + 1 


2 



Napomena: Slika prikazuje Sage radni list otvoren u WWW pregledniku s jednom izvršenom ćelijom 
i jednom praznom aktivnom ćelijom. U ovom dokumentu ćemo primjere koda davati onako kako bi 
izgledali u terminalskom sučelju gdje bi korisnikov upis slijedio nakon prompta sage :. Taj prompt ne 
smeta ni u ćelijama radnog lista, dakle primjeri iz ovog dokumenta se mogu slobodno kopirati u Sage 
grafičko sučelje sa ili bez prompta. Jedina razlika je da ukoliko više redaka iz ovog dokumenta kopiramo 
u jednu ćeliju, nećemo vidjeti medurezultate (ako ih trebamo, dodamo print ispred odgovarajućih 
redaka. Npr. kod iz ćelije sa slike će u ovom dokumentu izgledati ovako: 
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sage: 5-1 

4 

sage: 1+1 

2 


Otvaranje novih ćelija se izvodi stavljanjem miša između postojećih ćelija (ili na kraj radnog lista) i 
klikom u trenutku kad se pojavi plava horizontalna linija. Shift-klik će otvoriti tekstualnu ćeliju u 
koju stavljamo slobodni tekst i komentare. Brisanje ćelije se izvodi brisanjem sadržaja i onda pritiskom 
na backspace. 


Zadatak 1 

Kreirajte novu input ćeliju i izračunajte 2 + 3. Kreirajte zatim tekstualnu ćeliju s nekim tekstom 
iznad ove. Izbrišite obje ćeiije. 


Važno svojstvo Sage radnog lista je da svaku ćeliju možemo i naknadno editirati i onda ponovno izvršiti. 
(Editiranje tekstualnih ćelija iniciramo dvostrukim klikom miša negdje na ćeliji.) To je onda sve skupa 
sličnije radu u tabličnom kalkulatoru (spreadsheet) nego standardnom programiranju. Osim samog sa¬ 
držaja ćelija, trenutno stanje radnog lista (tzv. sesija, session ili proces) je određeno i vrijednostima 
varijabli u memoriji računala, a to općenito ovisi o redosljedu kojim su ćelije izvršene. Resetiranje 
tog stanja se izvodi putem izbornika Action pa Restart worksheet. (Sto je korisna operacija u 
trenucima kad se Sage radni list počne neočekivano ponašati.) Zbog toga je po završetku rada dobro 
prekontrolirati konzistentnost i reproducibilnost cijelog računa tako da resetiramo sesiju i izvršimo sve 
ćelije odjednom pomoću izbomikai Action pa Evaluate Ali. 

Svakom radnom listu pripada nezavisni Sage proces s nezavisnim varijablama. 


2.2 Elementarno računanje 

Sage se može koristiti kao obični kalkulator proizvoljne preciznosti: 

sage: 3*2+1 
7 


sage : sqrt( 9 ) 
3 


sage: factorial (22) 
1124000727777607680000 


sage: 13 "'2 3 

41753905413413116367045797 


sage: 13. "'2 3 

4.17539054134131e25 

Uočite različit tretman cijelih {integer) i realnih (floatingpoint) brojeva. Cijeli brojevi se uvijek tretiraju 
egzaktno, bez zaokruživanja ili odbacivanja nekih znamenaka. 
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Eksplicitno pisanje znaka množenja * se djelomično može izbjeći pozivanjem funkcije 
implicit_multiplication (True) , ali je i dalje nužno pisati * prije zagrada (kako bi 
bilo jasno da nije riječ o pozivanju funkcije). U ovom dokumentu se uvijek eksplicitno piše znak 
množenja. 

Za zapis velikih brojeva može se koristiti standardni Fortran/C zapis po kojem se npr. 3.2 • 10^ zapisuje 
ovako: 

sage: 3.2e4 

32000.0000000000 


Pridruživanje vrijednosti varijablama izvodi se znakom jednakosti “=”. 

sage : x = 4 
sage: print x 

4 


sage : x + 3 

7 

Standardne matematičke funkcije i konstante imaju uobičajena imena i ponašanje: 

sage: sin(pi) 

0 


sage: log(e) 

1 


sage: log(sin(pi)) 
-Infinity 


sage : oo + 1 

+Infinity 


Zadatak 2 

Izračunajte y/2\/e^. 


Primijetite da niste dobili numerički več simbolički rezultat. Sage če se uvijek, ako je moguće, odlučiti 
za egzatni simbolički rezultat. Ako nas ipak zanima numeričko približenje možemo koristiti funkciju 
numerical_approx (), koja se skraćeno zove i n {) : 

sage: n (sgrt (2*sqrt (e'^pi) ) ) 

3.10176639383605 

Funkcija n() ima i opcionalni argument digits kojim možemo zatražiti i veći broj decimala numerič¬ 
kog približenja: 

sage: numerical_approx(sqrt{2*sqrt(e^pi)), digits=50) 
3.1017663938360514951983183136487613918314124732866 

Funkcije često imaju brojne opcionalne argumente (keyword arguments) koji se navode proizvoljnim 
redom, ali svakako iza obaveznih argumenata. 


2.2. Elementarno računanje 
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Zadatak 3 

Odredite broj tt s točnošću od 500 deeimala: pridružite tu vrijednost varijabli x i ispišite je. Odre¬ 
dite sin (x). 

S kompleksnim brojevima radimo jednostavno. Samo treba imati na umu da je imaginarna jedinica 
reprezentirana velikim slovom I. 

sage: I ''2 

-1 


Zadatak 4 

Izračunajte -|- 1. 


Za upis matematičkih izraza u tekstualne ćelije stavljamo LaTeX kod unutar dolarskih znakova. Tako se 
upis 

$\alpha$ 

unutar retka prikazuje kao a. Veće jednadžbe koje trebaju stajati u posebnom redu upisuju se između 
parova dolarskih znakova u novom redu. Tako se: 

$$ E = \ganuna m c^2 $$ 

prikazuje kao 


E = 

Za formatiranje ispisa rezultata računa koristimo operator %, te standardne stringove za formatiranje 
(poznate iz sprintf() C funkcije): Npr. %5.3f formatira realni broj na 5 mjesta ukupno i 3 mjesta iza 
decimalne točke.) 

sage: print "%s = %10.8f ..." % ('Ludolfov broj', pi) 

Ludolfov broj = 3.14159265 ... 


2.3 Help sustav 

Da bismo pronašli potrebnu funkciju te način i primjere njene upotrebe služimo se slijedećim pristupima 
Sage dokumentaciji. Kao prvo, tu je TAB-nastavl janje (TAB-completion): započnemo li pisati ime 
neke funkcije, pritisak na tipku TAB dovršava pisanje njenog imena ako je nastavak jedinstven, a ako 
nije dobivamo popis svih mogućnosti (pa željenu odaberemo mišem ili kursorskim tipkama i tipkom 

Enter). 

Dokumentaciju konkretne funkcije dobijemo tako da nakon imena funkcije stavimo upitnik i onda pritis¬ 
nemo TAB (dobiveni tekst se može “odlijepiti” u poseban prozor pritiskom na pop-up link u njegovom 
desnom gornjem kutu). Iz dobivene dokumentacije možemo cut-and-paste-ati primjere u radni list. 
(Ukoliko umjesto jednog stavimo dva upitnika dobijemo cijeli ispis koda koji definira tu funkciju.) 
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Zadatak 5 

Otvorite novu ćeliju, utipkajte int pa pritisnite TAB. Odaberite funkciju integrate, dodajte 
upitnik ? pa ponovno stisnite TAB da dobijete osnovnu dokumentaciju. Prekopirajte jedan primjer 
upotrebe ove funkcije u istu ćeliju i izvršite je. 


Za pretraživanje po cijelom tekstu Sage dokumentacije postoji funkcija search_doc (), ali u praksi je 
efikasnije koristiti google koji relevantnije rezultate smješta na vrh. Kako je “sage” relativno česta riječ 
trik je prilikom pretraživanja koristiti kao ključnu riječ “sagemath” (što je alternativno ime tog softvera 
i nalazi se u WWW adresi Sage projekta). 


Zadatak 6 

Odredite numeričku vrijednost logaritma imaginarnog broja 16f u bazi 4, dakle log4(16f), tako da 
u dokumentaciji nađete kako se odabire drugačija baza funkcije log {) (defaultna baza je ona za 
prirodni logaritam tj. e = 2.718...). 


Zadatak 7 

Izvrijednite Eulerovu gama funkciju za z=l/2, dakle r(l/2), tako da pronađete u dokumentaciji 
odgovarajuču funkciju. 


Zadatak 8 


Proučite upotrebu funkcije sum 

) za zbrajanje matematičkih redova i izračunajte: 


(a) 1 + 3 + 5 + -- - + 61 


OD - 

X=1 


2.4 Poruke o greškama 

Ako se ogriješimo o matematička ili sintaktička pravila Sage če nam uzvratiti porukom o grešci. Klik 
mišem lijevo od vrha te poruke daje opširnije informacije (inače, drugi klik potpuno skriva poruku što 
se može koristiti i za skrivanje svih nepregledno dugačkih ispisa rezultata računa.) Za interpretaciju 
opširnije poruke potrebno je znanje Python programskog jezika, no ključna informacija je obično u 
zadnjem retku, koji je odmah vidljiv. 

sage: 1/0 

Traceback (most recent call last): 

ZeroDivisionError : rational division by zero 


sage: sin[2. 3] 

Traceback (most recent call last): 


2.4. Poruke o greškama 
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TypeError: 'Function_sin' object has no attribute '_getitem_' 

Početniku će te poruke izgledati nejasno, ali s vremenom će poprimati sve više smisla i treba ih uvijek 
čitati. Npr. gornja greška “object is unsubscriptable” je posljedica toga što smo za poziv funkcije 
sin umjesto okruglih zagrada upotrijebili uglate, koje služe za pristup elementima (tj. indeksima, 
subskriptima) polja i matrica. 
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POGLAVLJE 3 


Programiranje 


3.1 Liste i drugi spremnici 

3.1.1 Liste 

Liste su središnji objekti programiranja u Sageu i Pythonu. Srodne su poljima {array) iz standardnih 
programskih jezika (C, Fortran), ali imaju bitno više svojstava. Elementi lista mogu biti praktički bilo 
koji objekti. 

sage: ti = [1, 3, 5, 7, 9, 11]; ti 
[1, 3, 5, 7, 9, 11] 

sage: t2 = [sin, cos, tan]; t2 
[sin, cos, tan] 

sage: t3 = [[], ["jedan"], ti]; t3 
[[], ['jedan'], [1, 3, 5, 7, 9, 11]] 

Liste se mogu npr. zbrajati, množiti i mjeriti im se duljina: 

sage: print tl+t2 

[1, 3, 5, 7, 9, 11, sin, cos, tan] 

sage : 7 *[ 3 ] 

[3, 3, 3, 3, 3, 3, 3] 

sage: print len(tl) 

6 


Kako je Python objektno orijentirani računalni jezik, mnoge operacije se izvode putem tzv. metoda 
objekta, a to su objektu pridružene funkcije koje se pozivaju sintaksom ob j ekt. met oda (). I liste su 
objekti pa imaju niz metoda (vidi popis), od kojih je daleko najvažnija metoda append {) koja služi za 
dodavanje elemenata na kraj liste: 

sage: ti.append(13); ti 
[1, 3, 5, 7, 9, 11, 13] 

Pojedinim elementima liste se pristupa pomoču indeksiranja, gdje prvi element liste ima indeks 0, drugi 
ima indeks 1 itd. Pojedine elemente možemo promijeniti običnim pridruživanjem pomoču znaka jedna¬ 
kosti: 

sage: tl[l] = "tri" 

sage: print ti 

[1, 'tri', 5, 7, 9, 11, 13] 
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Ukoliko želimo pristupiti većem broju elemenata od koristi su tzv. rezovi liste (engl. sliče). Općenita 
sintaksa reza je lista [početak : kraj : korak] , gdje je početak uključen u rez, a kraj nije: 

sage : ti [2:6] 

[5, 7, 9, 11] 

sage: ti [2: 6 : 2 ] 

[5, 9] 

Ukoliko rez ide od samog početka ili sasvim do kraja liste, odgovarajuće indekse možemo izostaviti: 

sage : ti [2:] 

[5, 7, 9, 11, 13] 

Negativni indeksi nam omogućuju da brojimo od kraja liste ... 

sage: ti[-2:] # zadnja dva elementa 

[11, 13] 

...a negativan korak da rez ide unatrag tj. da se izvrne redosljed elemenata: 

sage: tl[::-l] # lista natraške 

[13, 11, 9, 7, 5, 'tri', 1] 


Zadatak 1 

Promijenite u gornjoj listi ti element "tri" (ne element s indeksom=3!) natrag u broj 3, ali ne 
eksplicitnom upotrebom indeksa i, već tako da “pronađete” indeks elementa "tri" ' 'pomoču 
metode ' 'index {). (Nije dozvoljeno upotrijebiti znamenku i). 


Za konstrukciju liste cijelih brojeva, vrlo je korisna funkcija range (početak, kraj, korak) 
gdje treba imati na umu da će konstruirana lista početi s elementom početak, ali će završiti s elemen¬ 
tom kraj-i. 

sage: range (3) # rezultira listom od tri elementa 

[ 0 , 1 , 2 ] 

Za konstrukciju liste realnih brojeva, vidi dolje odjeljak NumPy polja. 


3.1.2 Obuhvaćanje liste 

u praksi je kirurško mijenjanje pojedinih elemenata liste, kao u gornjem zadatku, vrlo rijetko. Listama se 
u pravilu rukuje kao cjelinama i najčešće se djeluje na sve njihove elemente. U standardnim računalnim 
jezicima u tu svrhu se koriste petlje, ali u Pythonu je najelegantnije koristiti tzv. obuhvaćanje liste (engl. 
list comprehension ): 

sage: [n+1 for n in ti] 

[2, 4, 6, 8, 10, 12, 14] 

Kombiniranjem obuhvaćanja liste s funkcijom range () možemo konstruirati najrazličitije liste. Npr. 
lista od pet slučajnih brojeva se može konstruirati ovako: 
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sage: [randomO for k in range(5)] 

[0.24117853317222127, 0.4504912135031752, 0.28942549153153674, 
0.86208567873778, 0.09328097538429081] 

Često se javlja potreba za izborom elemenata liste koji zadovoljavaju neki kriterij. To se može izvesti 
obuhvaćanjem liste uz dodatni uvijet: 

sage: t5 = range(l, 11) 

sage: [n for n in t5 if is_even(n)] # izbor parnih elemenata 
[2, 4, 6, 8, 10] 

Ovdje smo koristili funkciju is_even () koje je jedna iz velike porodice tzv. predikata. Predikat je 
termin iz matematičke logike koji označava funkciju koja poprima isključivo vrijednosti True ili False. 
Većina predikata definiranih u Sageu počinje s is_ (jer odgovaraju na pitanje “da lije Ispišimo ih 
(Koristimo prvo funkciju dir () koja ispisuje imena svih trenutno poznatih objekata i, drugo, to da su 
stringovi kao liste znakova pa možemo koristiti rezove na njima): 

sage: [p for p in dir() if p [ : 3] =='is_'' ] 

['is_2_adic_genus', 'is_32_bit', 'is_64_bit', 'is_AbelianGroup', 

'is_AbelianGroupElement', 'is_AbelianGroupMorphism', 

<.. . odrezan ispis . . .> 

'is_odd', 'is_optimal_id', ' is_pAdicField', 'is_pAdicRing', 

'is_package_installed', 'is_power_of_two', 'is_prime', 'is_prime_power', 

' is_pseudopriine', 'is_square', ' is_squarefree', ' is_triangular_nuinber'] 


Zadatak 2 

Konstruirajte listu svih prim-brojeva manjih od 100. Koliko ima prim-brojeva manjih od 
10^, 10^, 10®,...? Usporedite rezultate (i brzinu njegovog dobivanja) s ugrađenom funkcijom 
prime_pi(). 


Napomena: Mjerenje vremena potrebnog da se izvrši neka ćelija izvodi se stavljanjem %time u prvi 
red, a prekid računa koji traje predugo izvodi se putem izbornika Action pa Interrupt ili pritiskom 
na tipku Escape. 


%time 

sage: factor(2923003274661805836407421649242809468366377451741) 
1208925819614629174706189 * 2417851639229258349412369 
CPU time: 2.80 s, Wall time: 2.94 s 


Zadatak 3 

Izračunajte funkciju 7r(x) (broj prim-brojeva manjih od x) za x = 10^® statističkom metodom: 
Izaberite uzorak od n slučajnih cijelih brojeva između 1 i x i testirajte koliko ima prim-brojeva 
među njima. Iz tog udjela odredite 'k{x). Kolika veličina uzorka vam treba da bi relativna greška 
prema prime_pi(10^10) razumno često pala ispod 1%? 


3.1. Liste i drugi spremnici 
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Zadatak 4 

Kreirajte listu od 100 slučajnih brojeva iz intervala [ 0, 10), a zatim u toj listi ostavite 
samo brojeve koji se za manje od 0.02 razlikuju od cijelog broja. Naputak: Testirajte 
abs {x-round(x)). 


3.1.3 Ostali spremnici 

Osim lista, postoje i drugi spremnici (containers) za objekte. Kao prvo tu su tzv tuplovi (tuple) koji 
“izvana” izgledaju isto kao liste samo u okruglim zagradama. 

sage: tpl = (1, 3, 5, 7) 
sage: tpl [3] 

7 

Glavna razlika prema listama je da su tuplovi nepromjenjivi. Nije moguće niti promijeniti neki njihov 
element niti im dodati nove: 

sage: tpl[3] = 2 

Traceback (most recent call last): 

TypeError: 'tuple' object does not support item assignment 

Razlika u upotrebi između tuplova i lista je da su tuplovi obično heterogeni strukturirani skupovi u 
kojima pojedina mjesta u tuplu nose različita značenja, dok su liste obično homogeni skupovi istovrsnih 
objekata. Npr, koordinate neke točke je prirodno staviti u tupi (x, y), a ne u listu [x,y], ali niz točaka je 
prirodno staviti u listu takvih tuplova [(xl, yl), (x2, y2), ...]. 

Kad varijablama pridružujemo vrijednosti iz kraćih lista ili tuplova zgodno je koristiti tehniku raspaki¬ 
ravanja 

sage: x, y, z = (1, 2, 3) 

sage: print "Norma vektora = %f" % sqrt(x**2 + y**2 + z**2) 

Norma vektora = 3.741657 

Raspakiravanje se često koristi kod procesiranja listi čiji su elementi liste ili tuplovi. U slijedećem pri¬ 
mjeru pretvaramo listu 2D kartezijevih vektora u listu njihovih normi. Vidimo kako možemo kombinirati 
raspakiravanje po unutarnjoj s obuhvaćanjem 2D liste po vanjskoj dimenziji: 

sage: [sqrt(x**2 + y**2) for x,y in [(1,2), (3,4), (5,6), (7,8)]] 

[sqrt(5) , 5, sqrt(61), sqrt(113)] 
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Zadatak 5 

Koristeći raspakiravanje tuplova u kombinaciji s obuhvaćanjem liste pretvorite ovu listu vektora 

tt = [(ri, 6 »i), (r 2 , 6 ' 2 ), • • •] 

sage; tt = [{2.23, 1.11), (5.0, 0.93), (7.81, 0.88), (10.63, 0.85)] 

iz 2D polarnog sustava u kartezijev: 

[(xi,yi),(x2,y2),---] 

Nakon toga transformirajte nastalu listu kartezijevih vektora u listu kartezijevih vektora s cjelobroj¬ 
nim koordinatama, zaokruživanjem vrijednosti x i y koordinata pomoću funkcije round (). 


Daljnji spremnici koji nam stoje na raspolaganju su skupovi {set ), koji su skupovi različitih(!) objekata 
i s kojima možemo raditi standardne stvari poput unije, presjeka, komplementa ... 

sage: si = set ([10, 9, 10, 8, 7, 7 , 7 , 4]); si 

{4, 7, 8, 9, 10} 

(Uočite da se duplikati automatski izbacuju.) 

sage: sl.union(tl) 

(1, 3, 4, 5, 7, 8, 9, 10, 11, 13} 
sage: si.intersection (ti) 

{7, 9} 

sage: si.difterence (ti) 

{4, 8, 10} 

Zadnji, vrlo važan, spremnik kojeg ćemo spomenuti je rječnik (engl. dictionarj ). Riječ je o preslikava¬ 
nju key -> value, gdje je value bilo koji objekt, a key može biti tipično broj ili string (premda su i 
druge stvari dopustive kao key, npr tuplovi). 

sage: dl = {'a':l, 'c':3, 'b':2}; dl 
{'a' : 1, 'b' : 2, 'c' ; 3} 

sage: d2 = (1: sin, 2: cos, 3: tan}; d2 
{1: sin, 2: cos, 3: tan} 

sage: d3 = {'ime': 'pero', 'prezime': 'peric', 'spol': 'M'}; d3 
{'ime' : 'pero' , ' prezime' : ' peric' , 'spol' : 'M' } 

Redosljed elemenata u rječniku nema značenja. Pristup pojedinim elementima rječnika je moguć “in- 
deksiranjem” pomoću ključa: 

sage : dl[ 'c' ] 

3 

Na isti načinje moguče dodavati nove elemente u rječnik: 

sage : d2 [5] = log 

sage : d2.keys () 

[1, 2, 3, 5] 


3.1. Liste i drugi spremnici 
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Iteriranje po rječniku ne daje elemente već samo ključeve: 

sage: [it for it in d2] 

[1, 2, 3, 5] 


sage: [d2 [n] (1.) for n in d2.k:eys()] 

[0.841470984807897, 0.540302305868140, 1.55740772465490, 0.000000000000000] 


Zadatak 6 

Izračunajte (iteriranjem po rječniku, a ne ručno!) prosjek godina ljudi iz slijedećeg rječnika (rezul¬ 
tat treba biti decimalni broj): 

sage: ages = {'john' : 74, 'paul': 71, 'george': 71, 'ringo' : 73} 


Za kraj, i stringove možemo interpretirati kao liste znakova i tretirati ih kao takve 

sage: si = 'Samobor' 

sage : si [-3 :] 

' bor' 


3.1.4 NumPy polja 

Obične liste se rabe za razne namjene, uključujući simboličke račune, ali za numeriku su pogodnija tzv. 
NumPy polja. NumPy (Numerical Pjthon) je modul za python namjenjen baratanju s multidimenzional- 
nim brojčanim poljima kakva se često sreću u svim područjima znanosti. Riječ je o velikom modulu koji 
nije automatski prisutan u Pythonu ili Sageu, već ga treba učitati pomoću naredbe import 

sage: import numpy as np 

Ovdje se istovremeno s očitavanjem uvodi skraćeno ime np. Svi objekti i funkcije NumPy modula sada 
su dostupni kao np .<ime objekta>. 

Numpy polja su u računalu zapisana u neprekinutom slijedu memorijskih lokacija što omogućuje brži 
pristup no zbog toga svi elementi polja moraju biti istog tipa (integer, f loat, complex ...). 

sage: a = np.array (range ( 1 , 11 )); a # konverzija obične liste u NumPy polje 
array([ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]) 

Vidimo da je NumPy polje ispisano na drugačiji način nego obična lista, ali za razlikovanje tih dviju 
objekata izgled ispisa nije pouzdani kriterij. Neke operacije, poput print, ispisuju NumPy polje da 
izgleda kao obična lista 

sage:: print a 

[1 2 3 4 5 6 7 8 9 10] 

Ako postoji dvojba, možemo ispitati tip objekta komandom type () 

sage: type(a) 

<type 'numpy.ndarray' > 
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Pri konstrukciji NumPy polja izbor tipa objekata je automatski, ali moguće gaje i odrediti opcionalnim 
argumentom dtype: 

sage: b = np.array (range ( 1 , 11 ) , dtype=np.float64); print b 
[ 1. 2. 3. 4. 5. 6. 7. 8. 9. 10.] 

Tip objekata u NumPy polju je zapisan u dtype atributu polja *. 

sage: a.dtype 
dtype( ' int64') 
sage: b.dtype 
dtype (' float64') 

Numpy polja imaju dosta više metoda od običnih lista ... 

sage: dir(a)[-66:] 

['ali', 'any', 'argmax', 'argmin', 'argsort', 'astype', 'base', 'byteswap', 
'choose', 'clip', 'compress', 'conj', 'conjugate', 'copy', 'ctypes', 

'cumprod', 'cumsum', 'data', 'diagonal', 'dtype', 'dump', 'dumps', 'fill', 
'flags', 'flat', 'flatten', 'getfield', 'imag', 'item', 'itemset', 

'itemsize', 'max', 'mean', 'min', 'nbytes', 'ndim', 'newbyteorder', 

'nonzero', 'prod', 'ptp', 'put', 'ravel', 'real', 'repeat', 'reshape', 
'resize', 'round', 'searchsorted', ' setfield', ' setflags', 'shape', 'siže', 

'sort', 'sgueeze', ' std', 'strides', 'sum', 'swapaxes', 'take', 'tofile', 
'tolist', 'tostring', 'trače', 'transpose', 'var', 'view'] 

... no zbog efikasnosti implementacije nema upravo onih metoda koje imaju obične liste. To je zato jer 
je npr. umetanje elementa u listu vrlo “skupa” operacija koja uključuje brisanje i pisanje svih elemenata 
koji u memoriji dolaze nakon tog mjesta umetanja. Ako želimo raditi takve stvari upotreba NumPy polja 
nam ionako neće donijeti nikakve prednosti i treba koristiti obične liste. (Za konverziju NumPy liste u 
običnu koristimo funkciju list {).) 

Jedno od vrlo korisnih svojstava NumPy lista je da se većina matematičkih operacija prirodno distribuira 
po elementima liste (obične liste nemaju to svojstvo): 

sage: xs = np.linspace(1,10,3); xs 
array ([ 1. , 5.5, 10. ]) 

sage: xs + 1 

array ([ 2. , 6.5, 11.]) 

sage: sin(xs) 

array([ 0.84147098, -0.70554033, -0.54402111]) 

NumPy liste mogu biti i višedimenzionalne, no i one su u memoriji računala zapisane u jednodimen¬ 
zionalnom (ID) slijedu memorijskih adresa. Stoga je preoblikovanje elemenata ID NumPy u 2D polje 
proizvoljnog oblika računalno “jeftina” operacija, koja može biti od koristi 

sage: b=a.reshape(2,5); print b 
[ [ 1 2 3 4 5] 

[6 7 8 910]] 


sage: b.shape # "oblik" polja - broj redaka 1 stupaca 
(2, 5) 


* Atributi objekta su sve ono čemu se pristupa operatorom točkice .. Ranije spomenute metode su atributi koji se mogu 
pozivati kao funkcije. dtype je atribut koji je naprosto tip elementa polja. Za bolje razumijevanje svega ovog dobro je 
pročitati nešto o objektno-orijentiranom programiranju u Pythonu. 


3.1. Liste i drugi spremnici 
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Pristupanje pojedinom elementu višedimenzionalnog polja je moguće očitim indeksiranjem a[i][j]..., ali 
je efikasnije koristiti skraćeno indeksiranje: a[i, j,...]: 

sage: b[l,l] = 42. 

sage : print b 

[[ 1 2 3 4 5] 

[642 8 910]] 

Pažnja: radi efikasnosti, preoblikovanja i rezovi kroz NumPy polja ne kopiraju originalne elemente već 
samo manipuliraju pokazivačima na ista mjesta. Tako se npr. ova netom načinjena zamjena u listi b 
odražava i u listi a čijim preoblikovanjem je lista b nastala: 

sage : print a 

[ 1 2 3 4 5 6 42 8 9 10] 

Obične liste nemaju takvo ponašanje: 

sage: t5 = tl[2:4]; t5 # ti je obična lista 
[5, 7] 

sage: t5[0] = 'novi'; print t5 
['novi', 7] 

sage: ti 

[1, 3, 5, 7, 9, 11, 13] 

Inače, rezovi kroz dvodimenzionalne NumPy liste su elegantan način ekstrakcije redaka ili stupaca: 

sage: print b[:, 1 ] 

[2 42] 

Ukoliko trebamo liste realnih (float) brojeva, možemo koristiti NumPy inačicu funkcije range () (sin¬ 
taksa je ista kao za range ()) 

sage: np.arange (2 . , 9.9, 1.1) 

array([ 2. , 3.1, 4.2, 5.3, 6.4, 7.5, 8.6, 9.7]) 

Funkcija np . range ( ) je pogodna kad želimo specificirati korak liste. Kad želimo specificirati broj 
elemenata liste koristimo funkciju np.linspace(). 

sage: np.linspace {0, 10, 5) 

array([ 0. , 2.5, 5. , 7.5, 10. ]) 

Za daljnje detalje o NumPy poljima, vidi NumPy Tutorial. 
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Zadatak 7 

Za spajanje dvije jednako duge ID liste u 2D listu (zapravo listu tuplova) koristimo naredbu zip: 

sage; xs = [11, 12, 13, 14, 15] 

sage; ys = [-1, -2, -3, -4, -5] 

sage; points = zip(xs, ys); points 

[{11, -1), (12, -2), (13, -3), (14, -4), (15, -5)] 

Rastavite listu point s natrag na originalne dvije ID liste i to dvjema različitim metodama: 

1. Dva obuhvaćanja liste points (ova metoda nema veze s NumPy-em) 

2. Dva reza kroz listu points, koju prvo pretvorite u NumPy polje pomoću np . array. 


3.2 Kontrola toka izvršavanja 

3.2.1 Grananje if-then 

Sintaksa je 

if <condition>; 

<block 1> 
elif <condition>; 

<block 2> 
else; 

<block 3> 

Treba uočiti dvotočke prije blokova koda koji se uvjetno izvršavanju. Također je važno paziti na kon- 
zistentno uvlačenje blokova koda (tamo gdje bi u C-u bile vitičaste zagrade), jer je to jedina vodilja 
kompajleru da prepozna gdje su granice bloka. (Tradicionalno se uvlači za 4 mjesta. Editor u Sage 
ćelijama automatski radi to uvlačenje.) Npr: 

sage; if 2>3; 

....; print "veći je" 

....; elif 2==3; 

....; print "jednak je" 

....; else; 

....; print "manji je" 

manji je 

U uvjetima se mogu koristiti standardni logički operatori and, or, not, ..., a važno je uočiti da se test 
usporedbe jednakosti radi s dvostrukim znakom jednakosti ==. (Jednostruka jednakost = je rezervirana 
za operacije pridruživanja poput x=l.) Operatori uspoređivanja su standardni: ==, !=, <, <=, 

>, >=. 

sage: if 2>3 or 2<3: 

. . . . : print "nisu isti" 

nisu isti 


3.2.2 Petlje 

while petlja ima standardni oblik 


3.2. Kontrola toka Izvršavanja 


19 











Sage računalno okruženje za fizičare, Distribucija 1.1 


while <condition>: 
<body> 


gdje je condition logički uvjet. Npr. 

sage : k = 1 

sage: while k < 10: 

....: print k, 

....: k += 1 

123456789 

(Zarez nakon print naredbe sprečava automatski prelazak u novi red.) 

For petlja ima specifičnu sintaksu 

for X in <iterable>: 

<body> 


gdje je <iterable> lista, tuple, skup, rječnik, ... Za prijevremeno iskakanje iz petlje postoji komanda 
break. Npr. slijedeči algoritam pronalazi prvi cijeli broj čiji rastav sadrži više od pet različitih prostih 
faktora. 

sage: for i in range(l,le5): 

if len ( factor(i) )>5: 

. . . . : break 

sage : print i 

30030 


Alternativni i nešto elegantniji algoritam bi bio 

sage : i = 1 

sage: while len (factor(i)) <6 : 

....: i +=1 

sage: print i 

30030 


Zadatak 1 

Isprogramirajte petlju koja če pomnožiti sve brojeve od 1 do 7 (dakle, izračunati će 7!). Riješite 
zadatak jednom koristeći for petlju, a drugi put while petlju. 


3.3 Funkcije 

Mogućnost dehniranja novih funkcija je osnovna stepenica k naprednijem programiranju. Trebati ćemo 
razlikovati dvije vrste funkcija: 

1. Simboličke funkcije 

2. Python funkcije 

Ugrubo, simbolička funkcija dopušta više simboličkih manipulacija (poput integracije ih deriviranja). 
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ali ne može sadržavati kompleksne algoritme već samo jednostavne izraze Python funkcija može biti 
proizvoljno kompleksna, ali ne može se uvijek npr. simbolički integrirati. 


3.3.1 Simboličke funkcije 

Simboličke funkcije se definiraju na slijedeći prirodan način: 

sage : f{x) = x^2 

sage : f 

X I—> x^2 


S njima se mogu raditi sve uobičajene operacije: 

sage: f{ 2 ) 

4 

sage: f{{2*x+3)"2) 

(2*x + 3) ^^4 

sage: diff(f(x), x) # simboličko diferenciranje, vidi slijedeće poglavlje 
2*x 


sage: type (sin) 

<class 'sage.functions.trig.Function_sin'> 
sage: type (f) 

<type 'sage.symbolic.expression.Expression'> 

Naravno, funkcija može biti i funkcija od više variabli: 

sage: g{x, a) = x''a 

sage: g{ 4 , 2 ) 

16 


3.3.2 Python funkcije 

Python funkcije se definiraju korištenjem ključnih riječi def i return, te blokova koda koji su konzis- 
tentno uvučeni 

sage: def h (x) : 

. . . . : "Kvadriraj broj x . " 

. . . . : return x''2 


sage: print h(3) 
9 


Kao prvi red tijela funkcije može se, kao gore, staviti dokumentacijski string (tzv. docstring ) kojem 
se kasnije može pristupiti standardnim metodama pristupa dokumentaciji. (Dakle, h? če ispisati doku¬ 
mentacijski string.) 

Funkcija može imati i opcionalne argumente s defaultnom vrijednošču: 


^ Formalno, simboličke funkcije nisu po svom tipu “funkcije”, već “simbolički izrazi koji se mogu pozivati” (“callable 
symbolic expression’j'. Za detalje razlika između simboličkih i Python funkcija vidi ovdje) 


3.3. Funkcije 
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sage: def fun (x, n=l, b=0): 
. . . . : return x''n + b 

sage: fun (3, 4, 5) 

86 

sage: fun(3, b=5, n=4) 

86 

sage: fun(3) 


3 


(Uočite da kad smo eksplicitno imenovali argumente nismo morali paziti na njihov poredak.) 

Bilo što može biti argument funkeije. Najmoćnija stvar, obilato korištena u funkeionalnom pristupu 
programiranju, je da i same funkeije mogu biti argumenti funkeija: 

sage: def gun(f, x) : 

. . . . : "Komponiraj dvaput funkciju sa samom sobom" 

. . . . : return f{f(x)) 

sage: print gun(sin, 2) 

sin(sin (2)) 

sage: gun(log, 0.1) 

0.834032445247956 + 3.14159265358979*1 


Zadatak 1 

Isprogramirajte funkeiju fibProc (n) koja računa n-ti član Fibonaeeijevog niza (niz kod kojeg je 
svaki član definiran kao zbroj prethodna dva, a javlja se pri analizi idealizirane populaeije zečeva). 


Zadatak 2 

Isprogramirajte funkeiju f ibBinet ( n ) koja računa n-ti član Fibonaeeijevog niza putem Binetove 
formule 

— cos(n7r)<^“"' 

— 

gdje je Lp tzv. zlatni omjer (golden_ratio). Uvjerite se da dobivate dobre vrijednosti za neke 

n. 


3.4 Crtanje grafova 

3.4.1 2D grafovi 

Glavna naredba za ertanje 2D grafova je plot () 

sage: plot(sin(x), (x,0,2*pi), figsize=[ 4 , 2 ]) 
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Za kombiniranje plotova može se kao prvi argument staviti lista funkcija ... 

sage: plot([sin(x), log(x)], {x,0,2*pi), figsize=[4,2]) 



... ali je bolje naprosto zbrajati grafove pojedinih funkcija, što nam onda daje kontrolu nad svojstvima 
pojedinih linija (boja, debljina, ...). 

sage: P1 = plot(log(x), (x,0,2*pi)) 

sage: P2 = plot (sin (x), (x,0,2*pi), color='red', linestyle=''— 

....: thickness=3) 

sage: {P1+P2) . show {axes_labels= [ ' x' , ' f {x) ' ]/■ frame=True, axes=False, 

....: fontsize=16) 



X 


Premda to nije obavezno, gore je bilo logično svojstva pojedinih linija stavljati kao argumente funkcija 
plot {), a svojstva grafa kao cjeline u završni show (). Inače, korištenjem metode show () moguče 
je jednostavno ponovno iscrtavanje grafa uz promjenu nekih opcija: 

sage: P.show(frame=False, axes=True, aspect_ratio=2, fontsize=9) 


3.4. Crtanje grafova 
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«x) 



Funkcija plot () sama određuje raspon vrijednosti ordinate pogodan za crtanje zadanih funkcija. Me¬ 
đutim, ukoliko funkcija ima singularitet u području crtanja, to ispada loše: 

sage: plot ( 1 / (l-x"'2) , (x,-2,2)) 


2500 


2000 


1500 


1000 


500 


■1000 


Tada moramo eksplicitno ograničiti ordinatu opcijama ymin and ymax. (A može se i lijepo označiti 
položaj polova opcijom detect_poles.) 

sage: plot ( 1 / (l-x"'2) , (x,-2,2), detect_poles = 'shoinf' , yinin=-5, ymax=5) 
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Grafove često želimo upotrijebiti i izvan Sage radnog lista, npr. u nekom članku ili prezentaciji. Eks- 
portiranje grafova i drugih objekata postiže se uporabom funkcije save {). Format grafičke datoteke 
određenje ekstenzijom. Dopuštene ekstenzije su . png, .ps, .eps, .svg, and .sobj. 

sage: P.save (' /tmp/graf.png' ) 


Zadatak 1 

Pronađite ovu datoteku na disku i prikažite je pomoču nekog programa za pregledavanje slika. (Na 
Windowsima je potrebno pristupiti virtualnoj Linux mašini putem ssh protokola (vvinscp, putty,...) 
spajanjem na adresu na koju se spaja i Windows WWW browser, uz user=sage, pass=sage. Ako to 
nije moguče naprosto spremite sliku klikom desnim gumbom miša na nju, pa “Save Image As ...”.) 


Zadatak 2 

Nacrtajte graf funkcije log(x) između x=0.5 i x=1.5 bez ikakvih oznaka, okvira i osi, dakle samo 
liniju. 


Osim funkcija eksplicitno zadanih u obliku y = y{x) možemo crtati i funkcije zadane parametarski u 
obliku y = y{t), x = x{t). Npr: 

sage : var( 't' ) 

t 

sage: ff = (exp(cos(t)) - 2*cos(4*t) + sin (t/12)''5) 

sage: parametric_plot((ff*cos(t), ff*sin{t)), 

(t, 0, 2*pi) , fill=True, fillcolor='orange', figsize=[3,3]) 


3.4. Crtanje grafova 
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Za crtanje raznih potencijala i srodnih funkcija, zgodna je funkcija contour_plot () : 
sage : var( 'x y' ) 

sage: contour_plot {x''2-y'2, {x,-3,3), {y,-3,3), labels=True, 

label_colors='red',contours=[0,2,3,6], cmap='jet', colorbar=True) 

Aa 


T 



Zadatak 3 

Koristeći funkciju parametric_plot () i trigonometrijske funkcije nacrtajte kružnicu. 


Za crtanje podataka organiziranih u listu parova [[xi, ?/i], [x 2 ,y 2 ], ■ ■ ■] rabimo list_plot (): 

sage: data= [[0,0], [1,0.8], [2, 0.9], [3, 0.2], [4, -0.7]] 

sage: list_plot(data, pointsize=20, color='red', figsize=[ 4 , 2 ]) 
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Zadatak 4 

Uočite da se funkcija fibBinet (x) iz prošlog odjeljka može izvrijedniti i za vrijednosti ar¬ 
gumenta X koje nisu cjelobrojne. Nacrtajte graf te funkcije za -4 < x < 7 i superponi- 
rajte na njega (u drugoj boji) točke koje odgovaraju vrijednostima funkcije za pozitivne cjelo¬ 
brojne argumente [ (1, F_l) , (2, F_2) , ... (7, F_7) ]. Naputak: plot (...) 

+ list_plot (...). 


Zadatak 5 

Koristeći funkciju parametric_plot () i trigonometrijske funkcije nacrtajte kružnicu. 


3.4.2 3D grafovi 

Za crtanje 3D grafova stoje na raspolaganju dvije softverske biblioteke: 

• JMOL (default) traži da rade Java appleti u brovvseru (na Linuxu samo Sun Java funkcionira - 
mogući problemi u F-26). 

• Tachyon - raytracing paket, radi bez Jave, nije moguć interaktivni rad s crtežom (okretanje mišem) 

sage: y = var( ' y' ) 

sage: P2 = plot3d {x''4+y^4, {x, -2, 2), {y, -2, 2)); P2.show{) 

Klik mišem na sliku omogućuje mijenjanje smjera gledanja i zumiranje kotačićem miša. Desni klik 
otvara izbornik. Za “pravi 3D doživljaj” stavite dvobojne naočale pa onda desni klik -> Style -> Stere- 
ographic ->... 

Ukoliko nemate mogućnost prikazivanja Java appleta, koristite ‘tachyon’: 

sage: P2.show(viewer=' tacdvon' , figsize=[ 3 , 3] ) 
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sage: 
sage: 
sage: 
sage: 


L = plot3d (lambda x,y: 0, (-5,5), (-5,5), color="lightblue", opacity=0.8) 

P = plot 3d (lambda x,y: 4 - x^3 - y''2, (-2,2), (-2,2), color=' green' ) 

Q = plot3d (lambda x,y: x^3 + y^2 - 4, (-2,2), (-2,2), color=' orange' ) 

(L + P + Q).show(viewer=' tachyon' , figsize=[3,3]) 
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3.4.3 Matplotlib biblioteka 

Dosad rečeno je dovoljno za osnovno skiciranje. Funkcija plot {) i ostale navedene funkcije za 2D 
grafove su zapravo samo sučelje za tzv. matplotlib grafičku biblioteku. Za precizniju kontrolu nad 
crtanjem potrebno je izravno pozivati funkcije ove biblioteke. 

Matplotlib zapravo ne crta funkcije tj. linije več samo liste točaka (slično kao list_plot() gore), pa se 
crtanje linija dobiva iscrtavanjem i spajanjem gustih skupova točaka. Za kreiranje ovih lista točaka 
možemo koristiti NumPy funkcije linspace {) i logspace {). koristiti. 

sage: import matplotlib.pyplot as pit 
sage: import niampy as np 


Elementarni primjer 

sage: fig, ax = pit.subplots {) 
sage: xs = np.linspace (0, 2*np.pi) 
sage: ax.plot(xs, sin(xs)) 
sage: fig.savefig{ 'figl' ) 



Par komentara: 

• Za komunikaciju s Matplotlib bibliotekom koristimo modul pyplot (kojeg smo preimenovali u 
pit) koji omogučuje nešto pristupačnije sučelje 

• Kao granicu koristimo numeričku vrijednost za vr iz NumPy biblioteke. 

• Funkcija plot. subplots stvara dva objekta: sliku (figure) i panel (axis). Složene slike mogu 
imati više panela, vidi primjer dolje. 

• U čistom Pythonu, komanda za prikazivanje slike bila bi f ig . show (), ali u Sageu je za prikaz 
potrebno kreirati datoteku sa f ig. savef ig () koju onda Sage automatski prikaže u radnom 
listu 

^ Postoji i modul pylab sa sučeljem sasvim sličnim komercijalnom programu Matlab. Vidi usporedbu raznih Matplotlib 
sučelja. 

Primjere crtanja iz originalne matplotlib dokumentacije moguče je izravno kopirati u Sage, do na 
tu zamjenu show () sa savef ig ('fig') • Usput rečeno, savef ig () sprema datoteku sa slikom u 
\$HOME/ .sage/sage_notebook.sagenb/home/<username>/<kk>/cells/<nn>/imeslike.png, gdje je 
<kk> broj radnog lista vidljiv u URL-u (“WWW” adresi vidljivoj u WWW pregledniku), a <nn> je broj ćelije koji je vidljiv 


3.4. Crtanje grafova 
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• Najvažnije: sintaksa Matplolibove funkcije plot traži kao prve argumente posebno listu x- 
koordinata, i posebno listu y-koordinata: 

ax.plot([xl, x2, ...], [yl, y 2, ...], <opcionalni argumenti>) 

To je različito od srodne Sageove funkcije list_plot koja traži kao jedini argument listu parova 
koordinata: 

list_plot([[xl, yl], [x2, y2], ...], <opcionalni argumenti>) 

Za pretvorbu jedne vrste liste u druge, vidi zadatak na kraju odjeljka o NumPy listama. 

sage: fig, (axl, ax2) = pit.subplots(1, 2, sharey=True, figsize=[4, 3] ) 
sage: axl.plot(xs, sin(xs)) 
sage: axl.set_title (' sinus' ) 

sage: ax2.plot(xs, cos(xs), color='red', linestyle='— ') 
sage: ax2.set_title( 'kosinus' ) 
sage: fig.savefig{ ' figl' ) 


sinus 
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1 —I—I—I—I—r 
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i f 
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sage: 
sage: 
sage: 


sage: 


sage: 
sage: 
sage: 
sage: 
sage: 
sage: 
sage: 


xs = np . logspace(-2 . 0, 0.8, 100) # granice su log_10 (x) 

fig, ax = pit . subplots{figsize=[ 7 , 4 ]) # specifikacija dimenzija 

ax.plot(xs, sin(xs), color=' red' , linestyle='—', 

label='$\sin(x)$') # LaTeX oznake 

ax.plot(xs, sqrt(xs), 

'b-', label='$\sqrt{X } $' ) # skraćene oznake za boju i tip linije 
ax . set_xscale( ' log' ) 

ax . axhline( 0 , color='g', linewidth=l) # horizontalna linija na y=0 
ax . set_xlabel( 'x' , fontsize=14) 
ax . set_ylabel( 'f(x) ' , fontsize=14) 
ax . legend (loc="upper left") 

fig . tight_layout() # inače se ne vide cijele oznake na osima 

fig . savefig{ 'test' ) 


samo u tekstualnoj varijanti radnog lista (tipke "Text" ili "Edit" u notebooku), ali može se saznati i pomoću Unix shell 
komande f ind . -name " ime . png" iz direktorija ... <kk>/cells. 
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Zadatak 6 

Nacrtajte (Matplotlibom) 100 slučajno raspoređenih točaka (hint: np . random. rand () , te 
ax .plot ( . . ., linestyle=' None' , marker=' o' )) gdje su x i y koordinate tih točaka 
u intervalu (-1,1), a nacrtane su na dijagramu s x i y osima koje se protežu u intervalu (-2, 2). (hint: 
ax . set_xlim ( ) i ax . set_ylim ( )). 


Zadatak 7 

Odredite (x,y) koordinate vrha lijevog “krila” gornjeg “leptira”, dakle točke, negdje oko (x=3, y=3) 
koja je najudaljenija od ishodišta. 


Zadatak 8 

Nacrtajte kružnicu u kompleksnoj ravnini zadanu kompleksnim brojevima 

{z = pe*^(l- 1 I (/> G [0,27r),p = 1.3,r/= 0.2} 
te krivulju koja se dobije kad se ova kružnica podvrgne transformaciji Žukovskog: 


z ^ w 


1 

2 
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POGLAVLJE 4 


Matematika 


4.1 Simbolički izrazi 


Moć paketa za simboličku matematiku, poput Sagea, leži u sposobnosti manipulacije simboličkim izra¬ 
zima. Kao prvo, potrebno je na slijedeći način deklarirati varijable koje namjeravamo koristiti u simbo¬ 
ličkim izrazima ^ 

sage: var ('a b c x y z t') 

(a, b, c, X, y, z, t) 

Pomoću ovih varijabli sad izgrađujemo simboličke izraze: 

sage: ii = (b+a)^3; ii 
(a + b) "-3 

Sage ne provodi skoro nikakve operacije na izrazima dok to eksplicitno ne zatražimo. Recimo da želimo 
razviti gornji izraz koristeći binomni teorem. Za to služi funkcija expand ( ) 

sage: expand(il) 

a^3 + 3*a''2*b + 3*a*b''2 + b''3 

Funkcija expand (), kao i mnoge druge, se može alternativno upotrijebiti i kao metoda simboličkog 
izraza. 

sage: ii.expand() 

a^3 + 3*a''2*b + 3*a*b''2 + b''3 

U prvom pristupu expand doživljavamo kao funkciju ili operaciju, dok je izraz ii njen argument 
odnosno operand. To je način razmišljanja svojstven standardnom proceduralnom ili pak tzv. funkci¬ 
onalnom programiranju. 

U drugom pristupu izraz i 1 treba pak doživljavati kao objekt , u smislu tzv. objektno-orijentiranog (00) 
programiranja, aexpand{) je tzv. metoda što je naziv za funkciju koja je pridružena tipu objekta na 
koji djeluje 

* Eksplicitno deklariranje simboličkih varijabli je moguće izbjeći pozivanjem funkcije automatic_names (True), ali 
to nećemo koristiti. 

^ To onda omogućuje da metode istog imena rade različite stvari s različitim objektima (tzv. polimorfizam ). Kako je python 
OO jezik, takva sintaksa se obilato koristi u Sage-u i brojne funkcije se ni ne mogu koristiti na prvi način. Jedna od prednosti 
takvog pristupa je da elegantno možemo saznati popis svih funkcija koje rade nešto smisleno sa zadanim objektom, i to tako da 
nakon što stavimo točkicu poslije objekta stisnemo TAB tipku. Dobit ćemo popis svih metoda tog objekta. Ovo međutim 
ne funkcionira s netom upisanim izrazom u trenutnoj ćeliji, već samo s ranije definiranim izrazima (objektima) kojima smo 
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Da bi saznali što pojedina metoda radi, upišemo je nakon odgovarajućeg objekta i operatora toč¬ 
kice dodamo upitnik i stisnemo TAB. Pri upotrebi metode ne smije se zaboraviti na za¬ 
grade, koje su često prazne, ali nekad sadrže opcionalne argumente kojima modificiramo po¬ 
našanje metode. Ukoliko zaboravimo zagrade Sage ne poziva funkciju več samo ispisuje 
njeno puno ime poput "<metoda expand pridružena objektu 'Expression' koji 
je pohranjen na toj i toj adresi>". 

sage : ii. expand 

<built-in method expand of sage.syinbolic.expression.Expression object 

at 0x4e88200> 


Tek zagrade daju zahtjev interpreteru da dotičnu metodu i pozove tj. izvrši. 

Naravno ovakve jednostavne izraze možemo razviti i na ruke, dok računalo blista kad radi s velikim 
izrazima (sve dok stanu u memoriju računala) 

sage: 12 = (a +2*b + 3*c) ^3 * (x+y) ^3 
sage: i2.expand() 

a^3*x^3 + 6*a''2*b*x^3 + 12*a*b''2*x^3 + 8*b^3*x^3 + 9*a^2*c*x"'3 + 
36*a*b*c*x^3 + 36*b^2*c*x^3 + 27*a*c^2*x"'3 + 54*b*c^2*x''3 + 27*c^3*x^3 + 
3*a^3*x"'2*y + 18*a''2*b*x''2*y + 36*a*b^2*x''2*y + 24*b''3*x^2*y + 
27*a^2*c*x^2*y + 108*a*b*c*x^2*y + 108*b"'2*c*x"'2*y + 81*a*c^2*x''2*y + 
162*b*c''2*x^2*y + 81*c^3*x^2*y + 3*a^3*x*y''2 + 18*a''2*b*x*y^2 + 
36*a*b^2*x*y^2 + 24*b''3*x*y''2 + 27*a^2*c*x*y^2 + 108*a*b*c*x*y^2 + 
108*b^2*c*x*y^2 + 81*a*c^2*x*y''2 + 162*b*c''2*x*y^2 + 81*c^3*x*y''2 + 
a^3*y''3 + 6*a^2*b*y''3 + 12*a*b''2*y^3 + 8*b''3*y^3 + 9*a^2*c*y''3 + 
36*a*b*c*y''3 + 36*b''2*c*y''3 + 27*a*c^2*y''3 + 54*b*c''2*y^3 + 27*c''3*y^3 

Potenciranjem i razvijanjem gornjeg izraza dobivamo izraz od 550 članova ... 

sage: 13 = expand{i2^3) 
sage: len(i3) 

550 

... kojeg Sage s lakoćom faktorizira: 

sage: factor(i3) 

(a + 2*b + 3*c)^9*{x + y)^9 

Cesto je korisno izraz organizirati kao polinom u nekoj varijabli. Za to služi funkcija collect {) : 

sage: i4=i2.expand{).collect{y) 

sage: 14 

a^3*x^3 + 6*a''2*b*x''3 + 12*a*b''2*x^3 + 8*b^3*x''3 + 9*a^2*c*x"'3 + 
36*a*b*c*x^3 + 36*b^2*c*x^3 + 27*a*c^2*x"'3 + 54*b*c^2*x''3 + 27*c^3*x^3 + 

(a''3 + 6*a^2*b + 12*a*b^2 + 8*b^3 + 9*a^2*c + 36*a*b*c + 36*b^2*c + 

27*a*c^2 + 54*b*c^2 + 27*c^3)*y^3 + 3* {a"'3*x + 6*a^2*b*x + 12*a*b^2*x + 
8*b^3*x + 9*a^2*c*x + 36*a*b*c*x + 36*b^2*c*x + 27*a*c^2*x + 54*b*c^2*x 
+ 27*c^3*x) *y^2 + 3* (a^3*x^2 + 6*a^2*b*x''2 + 12*a*b''2*x^2 + 8*b''3*x^2 + 
9*a^2*c*x''2 + 36*a*b*c*x^2 + 36*b^2*c*x^2 + 27*a*c^2*x^2 + 54*b*c^2*x''2 
+ 27*c^3*x''2) *y 


pridjelili ime. Pridjeljivanje imena objektima se izvodi znakom jednakosti i korisno je ne samo zbog navedenog razloga već i 
inače radi lakšeg baratanja izrazima i kasnijeg referiranja na iste. 
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sage: len(i4) 

13 

(Uočite da collect {) ne pojednostavljuje koeficijente ^.) Za dobiti koeficijent uz neku potenciju 
neke varijable koristi se funkcija coef f icient () . Npr, koeficijent uz a® jest 

sage: 13.coefficient (a, 9) 

x^9 + 9*x''8*y + 36*x''7*y^2 + 84*x^6*y^3 + 126*x^5*y''4 + 126*x^4*y^5 + 
84*x^3*y^6 + 36*x^2*y''7 + 9*x*y^8 + y"'9 

Najsveobuhvatnija funkcija za pojednostavljivanje simboličkih izražaje simplif y_f uli () : 

sage: 15 = a/(1-a) + a/(l+a); i5 
a/(a + 1) - a/(a - 1) 


sage: 15.simplify_full() 

-2*a/ (a''2 - 1) 

Funkcija simplify_fuli () je kompozicija elementarnijih funkcija za pojednostavljivanje izraza. 
Jedna od tih elementarnijih funkcija je simplif y_trig () koja pri pojednostavljivanju rabi samo 
trigonometrijske identitete. 

sage : {sin(x)''4 + 2*sin(x)^2*cos(x)^2 + cos(x)^4). simplif y_trig () 

1 

Uočite da Sage neče naivno “pojednostaviti” izraze koji uključuju multifunkcije , kao na slijedečem 
primjeru, u kojem upoznajemo i važnu metodu subs {) koja služi za uvrštavanje vrijednosti varijabli i 
druge supstitucije u izrazima 

sage: 16 = log(a) + log(b) 
sage: i6.simplify_full() 
log(a) + log(b) 
sage: 16.subs(a=-l,b=-l) 

2*1*pi 

sage: 17 = log(a*b) 
sage: 17.subs(a=-l,b=-l) 

0 


.. end of output 


Zadatak 1 

Uzmite izraz (a + 6) ((c + yx)x + tx‘^) i algebarskim manipulacijama natjerajte Sage da ga prikaže 
u slijedećim ekvivalentnim oblicima (gdje redosljed faktora odnosno članova nije bitan): 

1. {a + h){tx + xy + c)x 

2. atx‘^ + ax‘^y + btx‘^ + bx‘^y + acx + bcx 


^ To se može postići na jedan od dva slijedeća zaobilazna načina 


4.1. Simbolički izrazi 


35 






Sage računalno okruženje za fizičare, Distribucija 1.1 


Zadatak 2 

Koristeći algebarske manipulacije pokažite da vrijedi 


sin^ X + cos^ X 

-= 1 — sin X cos X 

sin X + cos X 


Pazite na sintaksu: sin^ x se unosi kao sin(x)^3! 


Bilješke 

sage: sum (i4 . coef ficient (y, a) . f actor () *y''a for a in range(4)) 

(a + 2*b + 3*c) ^3*x''3 + 3* (a + 2*b + 3*c) ''3*x^2*y + 

3* (a + 2*b + 3*c) ^3*x*y^2 + (a + 2*b + 3*c) ^3*y''3 

sage: sum ( [it. factor () *y'eksp for it,ek;sp in i4 . coef f icients {y) ] ) 
(a + 2*b + 3*c) ^3*x^3 + 3* (a + 2*b + 3*c) ''3*x^2*y + 

3*(a + 2*b + 3*c)^3*x*y^2 + (a + 2*b + 3*c)^3*y^3 

sage: var('x y a b c' ) 

(x, y, a, b, c) 


4.2 Jednadžbe 

u Sageu kao znak jednakosti u jednadžbama stoji ==, jer je uobičajeni znak = rezerviran za pridjeljivanje 
vrijednosti simbolima odnosno pridjeljivanje imena izrazima. Rješavanje jednadžbi se izvodi funkcijom 

solve(): 

sage: sol = solve{2*x^2 - 1 == 0, x); sol 
[x == -l/2*sqrt(2), x == l/2*sqrt(2)] 

Valja primijetiti slijedeće: 

1. Prilikom poziva funkcije solve () treba eksplicitno naznačiti po kojoj varijabli se traži rješava¬ 
nje. 

2. Pronađena su oba rješenja kvadratne jednadžbe. 

3. Rezultat je ispisan u obliku liste jednadžbi. To omogućuje uvrštavanje rješenja u neki drugi izraz, 
ili u samu originalnu jednadžbu radi provjere. (Zato smo gore toj listi rješenja odmah pridjelili 
ime sol): 

sage: {x+y) .subs (sol [0]) 
y - l/2*sqrt(2) 


sage: (2*x^2 - 1 == 0 ).subs{sol[ 1] ) 

0 == 0 

Funkcija solve () može rješavati i sustave jednadžbi, ako se kao argumenti zadaju liste jednadžbi i 
odgovarajuća lista varijabli: 
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sage: sol2 = solve { [x"'2 + y"'2 == 1, x - 2*y == 0], [x, y] ) ; sol2 

[[x == -2/5*sqrt(5), y == -l/5*sqrt(5)], 

[x == 2/5*sqrt{5), y == l/5*sqrt(5)]] 

Da bi dobili numeričke vrijednosti ovih rješenja ne možemo jednostavno primijeniti funkciju 
numerical_approx () (tj. n ()) na ovu listu rješenja, jer n () nije metoda liste (a ni jednadžbe), 
već samo primitivnijih objekata: 

sage: n{sol2) 

Traceback (most recent call last): 

TypeError: unable to coerce to a ComplexNumber: <type 'list'> 

Stoga je potrebno nekako primijeniti n () samo na izraze s desne strane gornjih jednadžbi npr. 

sage : [ {x.subs(sol2[k][0]).n{), y.subs{sol2[k][1]).n()) for k in range { 2 )] 

[(-0.894427190999916, -0.447213595499958), (0.894427190999916, 

0.447213595499958)] 

Ovo isto možemo izvesti znatno elegantnije ^ koristeći raspakiravanje tupla, te metodu rhs () {right- 
hand side za pristup desnoj strani jednadžbe: 

sage : [ (x.rhs () . n (), y .rhs() . n ()) for x, y in sol2] 

[(-0.894427190999916, -0.447213595499958), (0.894427190999916, 

0.447213595499958)] 


Zadatak 1 

Riješite sustav jednadžbi 

= 1 , = 1 , 

Koliko ima rješenja? 

sage: var('x y' ) # nužno, jer smo ih gore u petlji redefinirali 

(x, y) 

sage: # edit here 


Ukoliko sustav jednadžbi ima beskonačno rješenja, dobit ćemo rješenje koje uključuje slobodni parame¬ 
tar ili više njih: 

sage: solve([x+y == 3, 2*x+2*y == 6],x,y) 

[[x == -rl +3, y == rl]] 

sage: solve([cos(x)*sin(x) == 1/2, x+y == 0],x,y) 

[[x == l/4*pi + pi*z38, y == -l/4*pi - pi*z38]] 

Gore je r<bro j> neki realni, a z<bro j> neki cijeli broj. Ovo radi samo sa sustavom jednadžbi, a ne 
i s jednom jednadžbom: 

Općenito, računalni kod koji uključuje duboko indeksiranje višedimenzionalnih lista, poput izraz [ 0 ] [ 1 ] je nepregle¬ 
dan i nepitonican. 

^ Naravno, postoji i Ihs () . 


4.2. Jednadžbe 
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sage: solve(sin(x)==0, x) 

[X == 0] 

Funkcija solve () daje simbolička (analitička) rješenja jednadžbi. Međutim, neka rješenja npr. jed¬ 
nadžbi viših stupnjeva nije moguče analitički zapisati. Npr, za slijedeču jednadžbu solve ( ) nam daje 
samo jedno trivijalno realno rješenje: 

sage: eq = 9*x^6 + 4*x''4 + 3*x^3 + x - 17 == 0 
sage: solve(eq, x) 

[x == 1, 0 == 9*x^5 + 9*x^4 + 13*x''3 + 16*x''2 + 16*x + 17] 

No znamo da ta jednadžba, šestog stupnja, mora imati šest kompleksnih rješenja. Ostalih pet se ne da 
zapisati drugačije nego kao numeričke (floating point) brojeve. Da bismo dobili ta rješenja koristimo 
metodu roots, gdje opcijom ring=CC tražimo rješenja u prstenu kompleksnih brojeva: 

sage: eq.roots (x, multiplicities=False, ring=CC) 

[-1.10301507262981, 1.00000000000000, 

-0.491102035999093 - 0.988331495372071*1, 

-0.491102035999093 + 0.988331495372071*1, 

0.542609572314000 - 1.05431152068711*1, 

0.542609572314000 + 1.05431152068711*1] 


sage: len (_) 
6 


Daljnji je problem da je i roots () zapravo analitički rješavač jednadžbi (koristi egzaktne, a ne nu¬ 
meričke metode), a neke jednadžbe se ne mogu analitički egzaktno riješiti, poput onih koje uključuju 
transcendentalne funkcije. 

sage: eq2 = 2 * arctan(x) == x^2 
sage: eq2.roots(ring=CC) 

Traceback (most recent call last): 

TypeError: Cannot evaluate symbolic expression to a numeric value. 

U tom slučaju moramo pribječi pravom numeričkom rješavanju, pomoču funkcije find_root (). 
Mala komplikacija s f ind_root () je da mu se treba dati interval u kojem traži rješenje i da če pronači 
samo jedno. Nakon toga moguče treba tražiti dalje u zadavanjem intervala koji isključuje več pronađena 
rješenja: 

sage: find_root(eq2, 0, 10) 

0.0 


sage: find_root(eq2, 0.1, 10) # isključujemo 0.0 

1.3717743420150883 

sage: find_root(eq2, 1.4, 10) # isključujemo i 1.37 

Traceback (most recent call last): 

RuntimeError : f appears to have no zero on the interval 


sage: find_root(eq2, -1., 1) 
-2.326813826671918e-21 
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Ovo zadnje rješenje je zapravo 0 jer funkcija f ind_root {) radi numeriku s konačnom preciznošću, 
koja po defaultu otprilike odgovara preciznosti double precision floating point varijabli u Fortranu ili 
C-u. Povećavanjem radne preciznosti pomoću opcionalnog argumenta xtol vidimo da se rješenje još 
više približava nuli: 

sage: find_root(eq2, -1., 1, xtol=le-30) 

4.3470345330696587e-32 

Sada se možemo uvjeriti da su dva pronađena rješenja zaista jedina, i to tako da skiciramo graf funkcije 

f{x) = 2 arctan(x) — 

(čije nultoćke su ekvivalentne rješenjima naše jednadžbe) i uočimo da siječe apscisu na samo dva mjesta. 

sage: plot (2*arctan (x) - x''2, (x, -3, 3), yinin=-4) 



Zadatak 2 


Riješite jednadžbu 

X 


tanx-= 0 . 


10 


Zadatak 3 

Pronađite numeričku vrijednost nekog kompleksnog rješenja jednadžbe sinx = 2 i provjerite uvr¬ 
štavanjem. 


Zadatak 4 

Riješite nejednadžbu — 12<0. 


4.2. Jednadžbe 
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Zadatak 5 

Pronađite pozicije lokalnog minimuma te lokalnog maksimuma gama funkcije r(x) koji su najbliži 
točki X = 0 


Zadatak 6 

Pronađite pozicije minimuma i maksimuma Besselove funkcije Ji(x) koji su najbliži točki x = 0. 


Bilješke 

sage: var ('a b c x y z t' ) 

(a, b, c, X, y, z, t) 

4.3 Matematička analiza 

Sage sadrži sve standardne operacije koje se uče u matematičkoj analizi, poput limesa, nizova, redova i 
diferencijalnog računa. 

4.3.1 Limesi 

Limesi se izvrijednjavaju funkcijom limit (): 

sage: limit(sin(x)/x, x=0) 

1 

sage: limit ( (l + l/x)''x, x=oo) 
e 


4.3.2 Razvoj u red 

Taylorov razvoj neke funkcije po nekoj varijabli oko neke točke do nekog reda radi funkcija taylor {) : 

sage: taylor(sin (x), x, 0, 5) 
l/120*x"'5 - l/6*x^3 t X 


Zadatak 1 

Odredite Taylorov razvoj funkcije f{x) = arctan(a:^ + 1) oko točke x = oo do šestog reda. 
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Zadatak 2 

Kolika je relativna pogreška koju radimo ako za izvrijednjavanje /(2) = arctan(5) koristimo red 
dobiven u gornjem zadatku? Da li bi greška bila manja da smo upotrebljavali razvoj do šestog reda, 
ali oko X = 0? 


Zadatak 3 

Za koji X greška razvoja funkcije f{x) = arctan(x^ + 1) do šestog reda oko točke x = 0 postaje 
jednaka greški istog razvoja, ali oko točke x = oo? 


4.3.3 Derivacije 

Deriviranje se radi funkcijom dif f () 

sage: diff(exp(2*x)*cos(3*x), x) 

2*cos(3*x)*e^(2*x) - 3*e^(2*x)*sin(3*x) 

Višestruko deriviranje: 

sage: diff(exp(2*x)*cos(3*x), x, 4) 

-119*cos {3*x) (2*x) + 120*e"'(2*x) *sin {3*x) 

Deriviranje izraza koji uključuje opču simboličku funkciju f (x) , korištenjem Leibnitzovog lančanog 
pravila: 

sage: f = function{ 'f' ) # simboličke funkcije treba deklarirati 

sage: diff(x^2*f(x), x) 
x^2*D[0](f)(x) + 2*x*f(x) 


4.3.4 Simboličko integriranje 

Simboličko neodređeno integriranje (integriramo izraz koji smo gore dobili deriviranjem): 

sage: integral (-3*6"'(2*x) *sin {3*x) + 2*e''(2*x) *cos {3*x) , x) . simplify_full () 
-4*cos {x) *e-'(2*x) *sin (x) ■'2 + cos (x) *e^ {2*x) 


Ovaj je izraz ekvivalentan onom gore, ali to ovdje nije eksplicitno. Ponekad je korisno eksperimentirati 
i s drugim algoritmima integriranja: 


sage : integral (-3*e'' (2*x) *sin{3*x) 


cos (3*x) *e'' (2*x) 


+ 2*e^(2*x)*cos{3*x), x, 

algorithm='sympy') 


Primijetite da se konstanta integracije podrazumijeva bez navođenja. 

Simboličko rješavanje određenih integrala: 

sage: integral ( x*log(x), x, 0, 1) 

-1/4 


4.3. Matematička analiza 
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Ukoliko integracija ovisi o nekom parametru, defaultni maxima algoritam za integriranje može tražiti da 
se izjasnimo o svojstivima tog parametrima o kojima ovisi rezultat integracije. Dodjeljivanje svojstava 
simboličkim varijablama izvodi se pomoću funkcije as sume (). Kako takve pretpostavke ne bi utjecale 
na kasnije računa treba ih što prije pobrisati s f orget {). 

sage: integral(y*sqrt(x*y), y, 0, l-x) 

Is x-l positive, negative, or zero? 

Traceback (most recent call last): 

TypeError: Computation failed since Maxima reguested additional constraints 
(try the command 'assume{x-l>0) ' before integral or limit 
evaluation, for example): 


sage: assume(x-l==0); integral(y*sqrt(x*y), y, 0, l-x).factor{); forgetO 
2/5*sqrt{-(x - l)*x)*(x - 1)^2 

U ovom konkretnom slučaju rezultat zapravo ne ovisi o predznaku od x, što maxima izgleda ne zna, ali 
sympy zna: 

sage: integral (y*sqrt (x*y) , y, 0, l-x, algorithm='' sympy' ) 

2/5*sqrt{x)* (-X + l)^(5/2) 


Zadatak 4 

Izračunajte neodređeni integral 


/ 


+ 3 

— X — 1 


dx 


i onda provjerite dobiveni rezultat deriviranjem i algebarskim manipulacijama. 


Zadatak 5 

Izračunajte (simbolički) dvostruki određeni integral 



Zadatak 6 

Izračunajte dvostruki neodređeni integral 



i provjerite rezultat deriviranjem i algebarskim manipulacijama. 


4.3.5 Numeričko integriranje 

Neki integrali su preteški ili se naprosto ne daju prikazati u zatvorenoj formi pa Sage vrača zadani izraz: 
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sage: integral(arctan(gamma{x)), x, 0, 1) 

integrale(arctan(gamma(x)), x, 0, 1) 

Tada nam ostaje numerička integracija pomoću numerical_integral () čiji rezultat je dan kao 
tupi. Prvi element tupla je rezultat integracije, a drugi procijenjena greška. 

sage: numerical_integral(arctan(gamma(x)) , 0, 1) 

(1.1020657425550975, 1.2235387620352894e-l4) 

Primijetite da se kod numeričke integracije ne navodi eksplicitno varijabla integracije već se ona odre¬ 
đuje automatski. Integrand ionako ne smije ovisiti nego o samo jednoj varijabli da bi se mogao numerički 
integrirati. 

Zadatak 7 

Izračunajte integral 



simbolički i numerički i usporedite rezultate. 


4.4 Linearna algebra 


u Sageu postoje vektori i matrice kao specijalni objekti, ali mi ćemo koristiti NumPy polja koja nude 
sličnu funkcionalnost zahvaljujući rutinama NumPy modula za linearnu algebru. (Za “domaću” Sage 
linearnu algebru vidi ovaj odjeljak.) 

sage: import numpy as np 
sage: import numpy.linalg as la 

sage: vi = np.array([1, 1, 2] ) 
sage: v2 = np . array((2,2, 4 ) ) 

Množenje skalarom je prirodno: 

sage: 3* vi 

array([3, 3, 6] ) 

Skalami i vektorski produkt vektora: 

sage: np.dot(vl, v2) 

12 

sage: np.cross(vl, v2) 
array ( [ 0, 0, 0]) 

Norma (“duljina”) vektora: 


4.4. Linearna algebra 
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sage: la.norin(v2) 

4.8989794855663558 


Množenje matrica te množenje matrice i vektora ide na prirodan način: 

sage: A = np. array { [ [ 1, 2, 1], [4, 3, 3], [9, 1, 7]]); A 


array 

( [ [1, 2, 1 

] , 


[4, 3, 3 

] , 


[9, 1, 7 

] ] ) 

sage: 

np. dot (A, 

A) 

array 

([[18, 9, 

14] , 


[43, 20, 

34] , 


[76, 28, 

61] ] ) 

sage: 

np. dot (A, 

vi) 

array 

([ 5, 13, 

24] ) 

Determinanta i inverz matrice: 

sage: 

la. d0t(A) 

# = 


-6.9999999999999982 


sage: la.inv(A) 
arraY([[-2.57142857, 
[ 0.14285714, 
[ 3.28571429, 


1.85714286, 
0.28571429, 
-2.42857143, 


-0.42857143], 
-0.14285714], 

0.71428571]]) 


Numerika po defaultu ide s double precision točnošću, pa kad provjeravamo da li je ^ ^ 
dobijemo egzaktno jediničnu matricu: 


sage: should_be_unit = np.dot{la.inv(A), A); 


array ( [[ 1.OOOOOOOOe+00, 

[ 2.77555756e-17, 

[ -1.22124533e-15, 


1.387778780-15, 
1.OOOOOOOOe+00, 
-7.771561170-16, 


should_b 0 _unit 
8.326672680-16], 
-2.775557560-17], 

1 . 000000000 + 00 ]]) 


Moguće je zatražiti od NumPyja da ispisuje ovakve male brojeve kao nule: 

sage: np.s 0 t_printoptions(suppr 0 ss=Tru 0 ) 


sage: should_ 

_b0_ 

_unit 


array([[ 1 ., 

0 . 

,, 0, 


[ 0., 

1. 

,, -0. 


[-0., 

-0 . 

1. 

.] ] 


= 1 ne 
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Zadatak 1 

Za datu matricu A definiramo svojstvene vektore (eigenvectors) v i njima pripadajuće svojstvene 
vrijednosti A (eigenvalues) kao rješenja matrične jednadžbe 

Av = Av . 

Pomoću funkcije la.eigO odredite svojstvene vrijednosti i svojstvene vektore matrice 

/ 2.3 4.5 \ 

V 6.7 -1.2 J ’ 

i provjerite da dobivena rješenja zaista zadovoljavaju gornju jednadžbu. 


Zadatak 2 

Kreirajte 3x3 matricu sa slučajnim realnim brojevima između 0 i 10. Invertirajte je i pomnožite s 
originalnom matricom te se uvjerite da dobijete jediničnu matricu. 


Dijagonalizacija matrice A je pronalaženje njenog rastava oblika 

A = PDP-^ 

gdje je D dijagonalna matrica. Takav rastav se također izvodi fukcijom la.eigO, koja daje svojstvene 
vrijednosti i svojstvene vektore matrice. Naime, stupci matrice P su upravo svojstveni vektori od A, a 
dijagonalna matrica D na dijagonali ima upravo odgovarajuće svojstvene vrijednosti: 

sage: A = np.array{ [ [3, 1], [1, 3]]); A 

array ( [ [3, 1], 

[1, 3] ] ) 

sage: evs, P = la.eig(A) 

sage: D = np.diag(evs); D 
array ( [ [ 4., 0 . ], 

[ 0., 2.]]) 

sage: np.dot(np.dot(la.inv(P), A), P) # provjeravamo P^-1 A P = D 
array([[ 4 ., 0 . ] , 

[0., 2. ]]) 

Neke matrice nije moguće dijagonalizirati, u slučaju čega će matrice P i D koje vraća la.eigO i 
dalje zadovoljavati 


AP = PD 


ali P neće biti invertibilna: 

sage: A = np.array{ [ [ 1, 1], [0, 1]]); A 

array ( [ [1, 1], 

[ 0 , 1 ]]) 

sage: evs, P = la.eig(A) 
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sage: np.dot(A, P) - np.dot{P, np.diag(evs)) 
array([[ 0 . , 0 . ] , 

[ 0 ., 0 . ]]) 

Neinvertibilnost od P se ogleda u ogromnim brojevima koje dobijemo kad “invertiramo” tu matricu: 

sage: print la.inv{P) 

[[ 1.OOOOOOOOe+00 4.50359963e+15] 

[ 0.OOOOOOOOe+00 4.50359963e+15]] 


4.5 Diferencijalne jednadžbe 


4.5.1 Simboličko rješavanje 


Osnovna Sage funkcija koja traži analitička rješenja običnih diferencijalnih jednadžbi je deso Ive (). 
Riješimo pomoću nje jednadžbu gušenog harmoničkog oscilatora 


dP 



+ 4y = 0 . 


Prilikom definiranja diferencijalne jednadžbe treba eksplicitno deklarirati simboličku nezavisnu vari¬ 
jablu (obično je to vrijeme t) i nepoznatu simboličku funkciju (ovdje je to y (t) ), pomoću funkcije 
function{) : 


sage: t = var( ' t' ) 
sage: y = function{ 'y' ,t) 

sage: desolve(diff{y, t, 2) + 2*diff(y, t) + 4*y == 0, y) 
(_K2*cos{sqrt (3) *t) + _Kl*sin(sqrt (3) *t) ) *e^ (-t) 


Napomene: 

1. Prisjetite se da je znak jednakosti u jednadžbama == 

2. Kao što je kod običnih jednadžbi trebalo eksplicirati po kojoj varijabli se traži rješenje, tako se 
kod diferencijalnih specificira po kojoj funkciji se fraži rješenje 

3. Rješenja uključuju nepoznate konsfanfe ( ki i k 2 ) koje freba odredili iz počelnih uvjefa 

Ukoliko želimo parlikularno rješenje za neke kon kr etne početne uvjete, zadamo te uvjete u dodatnom 
argumentu ics (= initial conditions), a koji je lista oblika ics = [ to, y{to), y'{to)]: 

sage: des = desolve(diff(y, t, 2) + 2*diff(y, t) + 4*y == 0, y, 

ics=[0, 1, 1]); des 

1/3*(2*sqrt(3)*sin(sqrt (3)*t) + 3*cos(sqrt(3)*t) )*e^ (-t) 

sage: plot(des, (0,5), figsize= [7,3]) + text('Guseni h.o.', (2,0.5)) 
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Zadatak 1 

Riješite simbolički diferencijalnu jednadžbu 


dt " 


uz početni uvjet y(0) = 1/2, provjerite rješenje uvrštavanjem, te ga nacrtajte za t u rasponu 0 < 
t < 1. 


4.5.2 Numeričko rješavanje 

Neke jednadžbe se ne mogu riješiti analitički, pa moramo pribječi numeričkom integriranju. Za to čemo 
koristiti funkciju odeint iz paketa scipy. integrate. 

sage: from scipy.integrate import odeint 
sage: import matplotlib.pyplot as pit 
sage: import numpy as np 


Kao prvi primjer integrirati čemo jednadžbu iz gornjeg zadatka. Za odeint tu jednadžbu treba tran¬ 
sformirati u oblik 


i korisnik treba definirati Python funkciju koja odgovara desnoj strani ove jednadžbe. 

U našem slučaju je f{y, t) = -|- 1 pa stoga imamo 

sage: def func(y, t) : 

....: return y**2 + 1 

Funkciju odeint treba pozvati s ovom funkcijom kao prvim argumentom, početnom vrijednošču y{0) 
kao drugim argumentom i listom vremenskih točaka za koje želimo odrediti položaj sustava kao trečim 
argumentom: 

sage: y0 = 0.5 

sage: ts = np.linspace(0, 1) 

sage: ts.shape 

(50, ) 

Funkcija odeint vrača listu položaja u traženim vremenskim točkama: 


4.5. Diferencijalne jednadžbe 
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sage: ys = odeint{func, y0, ts) 
sage: ys.shape 
(50, 1) 


sage: print "\n%8s %8s" % {'vrijeme', 'položaj') 

sage: print 18*'-' 

sage : for t,y in zip(ts[:5], ys[:5, 0]): 

....: print "%8f %8f" % (t, y) 

sage: print "%8s %8s'' % (' (...) ', ' (...) ') 

vrijeme položaj 


0.000000 
0.020408 
0.040816 
0.061224 
0.081633 
(...) 


0.500000 
0.525777 
0.552113 
0.579049 
0.606630 
(...) 


sage: fig, ax = pit.subplots{figsize=[4,3]) 
sage: ax.plot(ts, ys) 
sage: fig.savefig{ ' fig' ) 



Kao slijedeći primjer proučiti ćemo tzv. van der Polovu jednadžbu, koja je drugog reda: 


d?x 

dt'^ 


2\ dx 




Jednadžbe višeg reda se za numeričku analizu treba prirediti tako da se svaka pretvori u sustav od dvije 
jednadžbe prvog reda, što se izvodi uvođenjem nove funkcije y{t) = dx/dt uz pomoć koje gornju 
jednadžbu možemo pretvoriti u ekvivalentni sustav 


dx 
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Taj sustav sad treba zapisati u vektorskom obliku tako da poprimi strukturu 

^ = k = l,2 

gdje je dvokomponentni vektor y = {x,y) = [y [ 0 ] , y [ 1 ] ]. Sada funkcija koju priređujemo za 
odeint opet odgovara desnoj strani ove jednadžbe i mora kao rezultat vratiti dvokomponentni vektor 
/ = (/l,/2) 

sage: def func(y, t) : 

return [y[l], (1 - y[0]**2)*y [ 1] - y[0]] 

Početni uvjet je isto vektor (x(0), y(0)), a i rezultat će biti NumPy polje oblika (broj vremenskih 
točaka) x (broj nepoznatih funkcija). 

sage: ts = np.linspace (0,15, 150 ) 

sage: y0 = [1 . , 1. ] 

sage: ys = odeint(func, y0, ts) 

sage: ys.shape 

(150, 2) 


sage: fig, ax = pit.subplots(figsize=[7,3]) 
sage: ax.plot(ts, ys[:,0], label=' x(t) ' ) 
sage: ax.plot(ts, ys[:,l], label=' y(t)' ) 
sage: ax. legendo 
sage: fig.savefig( ' fig' ) 



Zadatak 2 

Nacrtajte graf brzine gornjeg rješenja Van der Polove jednadžbe y{t) = dx{t)/dt, ali ne u ovisnosti 
o vremenu t, već o položaju x{t) (tzv. fazni dijagram). 


Zadatak 3 

Riješite numerički diferencijalnu jednadžbu gušenog harmonićkog oscilatora. Nacrtajte funkcije 
y ( t ) i X ( t ), te, na posebnom grafu, fazni dijagram y (x). 


4.5. Diferencijalne jednadžbe 
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4.6 Statistika 


Sage ima implementirane osnovne statističke funkcije, ali nama će biti zgodnije koristiti paket 
scipy. stats čije mogućnosti su trenutno veće. 

sage: import scipy 
sage: import scipy.stats 
sage: import matplotlib.pyplot as pit 
sage: import numpy as np 


Kreirajmo prvo NumPy listu brojeva. NumPy liste imaju metode za izračun osnovnih statističkih veli¬ 
čina poput srednje vrijednosti, varijance, standardne devijacije itd. 

sage: xs=np.linspace ( 1 , 9,21); xs 

array([ 1. , 1.4, 1.8, 2.2, 2.6, 3. , 3.4, 3.8, 4.2, 4.6, 5. , 

5.4, 5.8, 6.2, 6.6, 7. , 7.4, 7.8, 8.2, 8.6, 9. ]) 


sage: print "%.2f H— %.2f" % {xs.mean(), xs.std()) 

5.00 +- 2.42 

Obične liste nemaju gornje statističke metode pa ukoliko nas zanima npr. srednja vrijednost brojeva u 
običnoj listi treba je prvo konvertirati u numpy listu funkcijom np . array {) ili možemo primijeniti 
Sage funkciju mean ( ) (dakle ne metodu). 

sage: mean([3, 5, 7, 9, 11]) 

7 

scipy.stats paket ima implementirane brojne statističke raspodjele: 

• norm - normalna (Gaussova) raspodjela 

• poisson - Poissonova raspodjela 

• gamma - gamma raspodjela 

• chi2 - raspodjela 

• t - studentova t raspodjela 


Svaka od tih raspodjela ima svoje parametre (npr. za Gaussovu raspodjelu su to srednja vrijednost i 
standardna devijacija). Oni se najčešće specificiraju putem opcionalnih argumenata, a točnu sintaksu 
saznajemo iz standardne dokumentacije. 

Najvažnije metode distribucija su: 

• pdf () - sama distribucija {probability distribution function) f{x) 

• cdf () - integral distribucije {cumulative distribution function) <b(x) = f{y)dy 

• rvs () - slučajni brojevi (rondom variates) koji slijede distribuciju 

Tako je npr. vrijednost u x = 0 normalne distribucije sa srednjom vrijednosti /r = 5 i standardnom 
devijacijom a = 1.5, f{x = 0, /i = 5, cr = 1.5) dan s: 

sage: scipy.stats .norm.pdf (0, loc=5., scale=1.5) 

0.0010281859975274036 
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Skica distribucije i njenog integrala: 

sage: xs = np.linspace (-2, 12 ) 

sage: fig, ax = pit.subplots{figsize=[5,3]) 

sage: ax.plot(xs, scipy.stats.norm.pdf{xs, loc=5 . , scale=1.5), label='pdf') 
sage: ax.plot(xs, scipy.stats.norm.cdf{xs, loc=5 . , scale=1.5), 'r— 

.. . . : label='cdf') 

sage: ax.legend(loc=' upper left' ).draw_frame( 0 ) 
sage: fig.savefig{ 'fig' ) 



Izračunajmo vjerojatnost da se slučajna varijabla nade unutar jedne standardne devijacije a od srednje 
vrijednosti Ona je dana integralom distribucije od fj, — a do ^ + a. Kako na raspolaganju imamo 
funkciju cdf () koja daje integral od — oo do x, treba samo oduzeti dva takva integrala (koristimo i to 
da su defaultne vrijednosti loc=0 i scale=l): 

sage: scipy.stats .norm. cdf(1) - scipy.stats. norm .cdf(-1) 

0.68268949213708585 


(To je čuvenih 68% vjerojatnosti.) 


Zadatak 1 

Slučajna varijabla X slijedi normalnu raspodjelu sa srednjom vrijednosti /i = 50 i standardnom 
devijacijom a = 2. Kolika je vjerojatnost da se X nade između 47 i 54? 


Sad čemo metodom rvs () izgenerirati uzorak od 10000 brojeva raspodjeljenih po normalnoj raspodjeli 
s// = 5icr = 1.5i testirati da su srednja vrijednost i standardna devijacija prema očekivanjima: 

sage: pts=scipy.stats.norm.rvs(loc=5 . , scale=1.5, size=le4) 
sage: print "broj točaka = %i'' % len(pts) 
broj točaka = 10000 

sage: print "srednja vrijednost = %.3f" % pts.meanO # rel tol 0.01 
srednja vrijednost = 4.998 

sage: print "standardna devijacija = %.3f" % pts.stdO # rel tol 0.03 
standardna devijacija = 1.493 

Za crtanje histograma ovog uzorka koristimo Matplotlibovu hist () funkciju, čiji opcionalni argument 
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bins kontrolira broj “binova” dakle broj intervala apscise u kojima se prebrojavaju točke: 

sage: fig, ax = pit.subplots(figsize=[3 . 5, 2] ) 
sage: ax.hist(pts, bins=20, normed=True) 
sage: fig.savefig{ 'fig' ) 



Ovo je bilo za jedan uzorak. Listu npr. standardnih devijacija za 5 nezavisnih uzoraka možemo konstru¬ 
irati ovako: 

sage: [scipy.stats.norm.rvs{loc=5, scale=1.5, size=le4).std{) for k in 

range(5)] # rel tol 0.05 

[1.503340173922008, 1.5040004762159596, 1.4785845284471362, 
1.4906248759178866, 1.5136584591557063] 


Zadatak 2 

Kolika je vjerojatnost da mjerenje slučajne varijable koja slijedi normalnu razdiobu sa srednjom 
vrijednošću ^ i standardnom devijacijom a odstupi od srednje vrijednosti za više od Scr? 

Važna raspodjela je tzv. x^-raspodjela 




(^2)./2-lg-xV2 

2'^/2r(i//2) 


Vidimo da raspodjela, pored argumenta x, ima samo jedan parametar: v = Af = broj stupnjeva slobode 
{degrees offreedom) 

Srednja vrijednost raspodjele jednaka je broju stupnjeva slobode: 

sage: scipy.stats.chi2.rvs(3, size=le5).mean() # rel tol 0.01 

3.0033578764955071 


Jedno od svojstava raspodjele je da integral repa te raspodjele (s jednim stupnjem) od 1, 4, 9 do oo 
daje vjerojatnost da mjerenje slučajne varijable raspodjeljene po Gaussovoj raspodjeli sa standardnom 
devijacijom a odstupi od srednje vrijednosti za cr, 2cr, 3cr,.... 

sage: [l-scipy.stats.chi2.cdf{eps^2, 1 ) for eps in [1,2,3]] 
[0.31731050786291404, 0.045500263896358528, 0.0026997960632602069] 


® Naziv tog parametra postaje jasan u kontekstima određivanja dobrote prilagodbe i statističkog testiranja hipoteze. 
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Zadatak 3 

Središnji granični teorem kaže da su srednje vrijednosti dovoljno velikih uzoraka neke raspodjele 
raspodjeljene prema normalnoj (Gaussovoj) raspodjeli bez obzira na to kakva je originalna raspo¬ 
djela iz koje te uzorke uzimamo. U to ćemo se uvjeriti na primjeru raspodjele. 

1. Nacrtajte tu raspodjelu za df=3 stupnja slobode i uvjerite se daje asimetrična oko srednje 
vrijednosti. 

2. Uzmite jedan uzorak od milijun točaka te raspodjele i uvjerite se daje standardna devijacija 
a = \/2 ■ df 

3. Uzmite 500 uzoraka od po 400 točaka svaki i napravite listu srednjih vrijednosti tih uzoraka. 
Uvjerite se daje standardna devijacija vrijednosti u listi (t/\/ 400, kako traži teorem. 

4. Nacrtajte histogram srednjih vrijednosti u listi i uvjerite se (superponiranjem krivulje nor¬ 
malne razdiobe) da su one zaista raspodjeljene prema normalnoj razdiobi srednje vrijednosti 
/i = d/ = 3 i standardne devijacije cr/VdOO 


Za ove iste primjere, ah izvedene pomoću Sage funkcija, vidi odjeljak Statistika - alternative. 

Bilješke 

4.7 Prilagodba funkcije podacima 


često je potrebno neke podatke dobivene npr. mjerenjima grahčki prikazati i približno opisati nekom 
funkcijom (tzv. prilagodba ih “htanje”)- 


sage: data 


[[ 0 . 2 , 


0.1], [0.35, 


0 . 2 ] , 


[ 1 , 0 . 6 ], 


[ 1 . 8 , 


0.9], [2.5, 1.3], 

[4, 2.06], [5, 2.6]] 


:: sage: hst_plot(data, hgsize=[3,3]) 

2.5 - 

2 : 

1.5 - 

1 - 

0.5 - 

_._i___I_._I___I___L 

1 2 3 4 5 

Klasična situacija je potreba za prilagodbom pravca f{x) = a + bx metodom najmanjih kvadrata (li¬ 
nearna regresija). Pogledajmo prvo kako bismo sami implementirah standardne formule za metodu 
najmanjih kvadrata: 


4.7. Prilagodba funkcije podacima 
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sage: xs = [x for x,y in data] 
sage: ys = [y for x,y in data] 
sage: xs2 = [x^2 for x,y in data] 
sage: ys2 = [y^2 for x,y in data] 
sage: xys = [x*y for x,y in data] 
sage: n = len (data) 

sage: delc = n*sum{xs2) - sum{xs)''2 

sage: b = (n*sum(xys) -sum (xs) *sum (ys))/delc 

sage: a = (sum(xs2) *sum (ys) -sum (xs) *sum {xys))/delc 

sage: sigsq = sum ([(y-b*x-a) ^2 for x,y in data])/(n-2) 

sage: erra = sqrt( sum (xs2)*sigsq/delc) 

sage: errb = sqrt(n*sigsq/delc) 


sage: 

print 

"a = %. 

.3f +- 

%.3f" % 

(a. 

erra) 

a = 0 

.020 +- 

0.023 





sage: 

print 

"b = %. 

.3f +- 

%.3f'' % 

(b. 

errb) 


b = 0.513 +- 0.008 

Sage sadrži funkciju find_fit () koja radi otprilike to isto, ali je općenitija u smislu da krivulja 
koju prilagodavamo može biti zadana bilo kojim matematičkim izrazom. Pogledajmo prvo prilagodbu 
pravca: 

sage : var( 'a b' ) 

(a, b) 

sage: model(x) = a + b*x 


sage: fit = find_fit(data, model, solution_dict=True); fit 
{b: 0.5130561168359806, a: 0.020159524186939004} 

Ovdje smo zatražili ispis rješenja u obliku rječnika kako bismo ga lako uvrstili u model putem metode 
subs {): 

sage: list_plot(data, figsize=[3,3]) + plot(model(x) .subs (fit), (x, 0, 5), 

. . . . : color='red') 
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Zadatak 1 

Donjim podacima iz liste data2 prilagodite polinom trećeg reda, te funkciju f{x) = asin(x) te 
nacrtajte sve najednom grafu. Radi lakšeg snalaženja neka krivulje budu različitih boja. Izračunajte 
zbroj kvadrata odstupanja Yliiyi ~ oba slučaja. 


sage: data2 = [[0, 0.1], [0.3, 0.4], [1, 0.7], [1.8, 1], [2.5, 0.5], 

-: [4, -0.6], [5, -1.0]] 

Nažalost, find_fit() ne daje greške parametara pravca a i 6. Za to trebamo koristiti naprednije funkcije za 
prilagodbu. Zgodna takva funkcija je leastsg () iz paketa scipy. optimize. Ona kao svoj prvi 
argument traži funkciju koju minimiziramo (kod metode najmanjih kvadrata to je odstupanje funkcije 
u točki od ordinate točke, yi — f{xi), ah ne kvadrirano jer leastsg () sam radi kvadriranje!). Tu 
funkciju treba definirati kao Python funkciju. Npr. za prilagodbu pravca funkcija koju minimiziramo 
treba vračati listu odstupanja pravca od točaka: 

sage: def distances (p, data): 

....: return [(y - p[0] - p[l]*x) for x,y in data] 

Prvi argument ove funkcije je lista p [ ] parametara koje prilagodavamo. Početne vrijednosti tih para¬ 
metara specificiraju se u drugom argumentu funkcije“leastsq()“. 

Treći argument* leastsq() je tuple s eventualnim dodatnim argumentima funkcije koju minimizi¬ 
ramo (prvi argument te funkcije je uvijek lista parametara koje prilagodavamo). Ovdje ćemo taj dodatni 
argument iskoristiti da proslijedimo toj funkciji same podatke (data) na koje prilagođujemo funkciju. 

Po defaultu, leastsg {) isto vrača samo tupi s listom prilagođenih vrijednosti parametara (bez gre¬ 
šaka) i zastavicom (flag) koja označava tip pronađenog rješenja ih eventualne probleme: 

sage: import scipy 

sage: scipy.optimize.leastsq(distances, [l.,l.], (data,)) 

(array{[ 0.02015952, 0.51305612]), 3) 

Ukoliko želimo i greške koristimo opcionalni argument full_output=True, koji nam daje i tzv. 
matricu kovarijanci (cov_x dolje) iz koje odredimo greške parametara na slijedeći način 

sage: p_final, cov_x, infodict, msg, ier = scipy.optimize.leastsg( 

....: distances, [l.,l.], (data,), full_output=True) 

sage: vr = scipy.dot(distances(p_final, data), distances( 

....: p_final, data))/(len(data)-len(p_final)) 

sage: p_errs = sqrt(scipy.diagonal(vr*cov_x)) 

sage: print "p[0] =%8.3f +- %.4f" % (p_final[0], p_errs[0]) 

p[0] = 0.020 +- 0.0231 

sage: print "p[l] =%8.3f +- %.4f" % (p_final[l], p_errs[l]) 
p[l] = 0.513 +- 0.0085 

leastsq () također omogućuje i prilagodbu složenijih funkcija. Npr. možemo se pitati trebamo li 
gornjoj funkciji dodati i kvadratni član: 

’ Gornji način određivanja greške parametara se može prihvatiti kao recept, no za one koji žele znati evo objašnjenja: 
Kad se minimizira ispravno normalizirana funkcija onda dijagonalni elementi matrice kovarijanci izravno daju greške pa¬ 
rametara. No, za običnu prilagodbu metodom najmanjih kvadrata, minimizira se kvadrat apsolutne greške pa je potrebno 
renormalizirati matricu kovarijanci varijancom mjerenja, a procjenitelj te varijance je upravo var=(suma kvadrata 
odstupanja) / (broj stupnjeva slobode). Inače, kad se tako vrijednost kvadrata odstupanja iskoristi za procjenu 
greške mjerenja ista se više ne može jednostavno iskoristiti za procjenu dobrote fita. (Vidi diskusiju u Bevingtonu ispod Eq. 
(6.15).) 
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sage: def distances (p, data): 

return [{y - p[0] - p[l]*x - p[2]*x^2) for x,y in data] 

sage: p_final, cov_x, infodict, msg, ier = scipy.optimize.leastsg( 

distances, (data,), full_output=True) 

sage: vr = scipy.dot(distances(p_final, data), distances(p_final, 

. . . . : data))/(len(data)-len{p_final)) 

sage: p_errs = sqrt(scipy.diagonal(vr*cov_x)) 
sage: for k in range ( len (p_final)): 

print "p[%i] =%8.3f +- %.4f" % (k, p_final[k], p_errs[k]) 
p[0] = 0.027 +- 0.0350 

p[l] = 0.502 +- 0.0373 

p[2] = 0.002 +- 0.0071 

Činjenica da je greška paramera p [ 2 ] znatno veća od vrijednosti samog parametra, sugerira da je 
kvadratni član suvišan. (Štoviše, vidimo daje i slobodni član p [ 0 ] od malog značaja.) 

Bilješke 
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POGLAVLJE 5 


Fizika 


5.1 Mehanika 


Zadatak M-1 

Neka na česticu u ravnini djeluje potencijal U{x,y) = —Nacrtajte U(x,y) i pomoću 
contour_plot () i pomoću plot3d {). Nadalje, rješavajući numerički Nevvtonovu diferenci¬ 
jalnu jednadžbu gibanja, pronađite putanju čestice, uz početne uvjete x(0) = y{0) = —1, ^^(O) = 
2, Vy{0) = —3, od trenutka f = 0 do trenutka t = 1.43. Nacrtajte tu putanju, a onda je pokušajte 
superponirati na oba gornja dijagrama potencijala. Igrajte se malo s početnim uvjetima i uvjerite 
se da je sve OK. Izvrijednite ukupnu energiju za razna vremena i uvjerite se da je sustav zaista 
konzervativan tj. daje energija konstantna. 


Zadatak M-2 


Formule za nerelativističku i relativističku kinetičku energiju tijela mase m koje se giba brzinom v 
su 


Knr = —mv^ 


Kj-ei = me? 


{ , ‘ 

y Y^i — u2/c2 



1. Definirajte odgovarajuće funkcije knr (v) i krel (v) te usporedite na istom dijagramu 
ponašanje tih funkcija za brzine od 0 do 3 • 10^ m/s za jediničnu vrijednost mase. Označite 
veličine i njihove jedinice na koordinatnim osima! 

2. Odredite brzinu pri kojoj je greška nerelativističke formule tisučinku promila. 

3. Ako ste definirali funkcije kao Python funkcije, ponovite zadatak sa simboličkim funkcijama 
i obratno. 

4. Da li se knr (v) i krel ( v) slažu za male vrijednosti brzine? Ako ne, korigirajte krel {) 
da postignete slaganje. 


5.1.1 Problem tri tijela 

Numeričkim rješavanjem Newtonovih jednadžbi simuliramo gibanje tri tijela koja gravitiraju, a čije 
je gibanje ograničeno na x-y ravninu. Radimo s jediničnim masama i G=1 vrijednošću Nevvtonove 
gravitacijske konstante. Za tri tijela u dvodimenzionalnoj ravnini imamo 6 diferencijalnih jednadžbi 
drugog reda koje pretvaramo u 12 diferncijalnih jednadžbi prvog reda. 
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sage 


def system(t, y): 

# y_i = {xl, yl, x2, y2, x3, y3 

#i= 012345 

# vxl, vyl, vx2, vy2, vx3 

# 6 7 8 9 10 

return [y[6], y[7], y[8], y[9], 
-((y[0] - y[2])/((y[0] - y[2])** 
(y[0] - y[4])/((y[0] - y[4])**2 
-((y[l] - y[3])/({y[0] - y[2])** 
(y[l] - y[5])/((y[0] - y[4])**2 
(y[0] - y[2])/((y[0] - y[2])**2 
(y[2] - y[4])/((y[2] - y[4])**2 
(y[l] - y[3])/((y[0] - y[2])**2 
(y[3] - y[5])/((y[2] - y[4])**2 
(y[0] - y[4])/((y[0] - y[4])**2 
(y[2] - y[4])/( (y[2] - y[4] )**2 
(y[l] - y[5])/((y[0] - y[4])**2 
(y[3] - y[5])/((y[2] - y[4])**2 

] 


vy3) 

11 

y[10], y[ll]. 


+ (y[l] 
(y[l] - 
+ (y[i] 
(y[l] - 
(y[l] - 
(y[3] - 
(y[l] - 
(y[3] - 
(y[i] - 
(y[3] - 
(y[i] - 
(y[3] - 


- y[3]) 
y[5])** 

- y[3]) 
y[5])** 
y[3])** 
y[5])** 
y[3])** 
y[5])** 
y[5])** 
y[5])** 
y[5])** 
y[5])** 


+*2)**1,5) 
2)**1.5, 
**2)**1,5) 
2 )** 1 . 5 , 

2)**1.5 - 
2 )** 1 . 5 , 

2)**1.5 - 
2)**1.5, 

2)**1.5 + 
2)**1.5, 

2)**1.5 + 

2)**1.5, 


Za početne uvjete uzimamo neke od 15 periodičkih orbita navedenih u Tablici 1 članka Šuvakov i Dmi- 
trašinović, Phys. Rev. Lett. 110, 114301 (2013), dostupno online: arXiv:1303.0181. 

sage: figure8=[[-1, 0, 1,0,0, 0, 0.347111,0.532 728, 0.347111,0.532728,-0.694 222, 
....: -1.06546],[-1.05783,1.05784,-0.522919,0.522917],6.32445] 

sage: goggles=[[-1,0,1,0,0,0, 0.0833,0.127889,0.0833,0.127889,-0.1666, 

....: -0.255778], [-1.04973,1.04972,-0.738248,0.738249],10.4668] 

sage: butterfly2=[[-1,0, 1, 0, 0, 0, 0.392 955, 0.097579, 0.392 955, 0.097579,-0.78591, 
....: -0.195158], [-1.07444,1.07449,-0.174197,0.17414],7.00391] 

sage : yinyang2a=[[-1, 0, 1,0, 0, 0, 0.416822,0.330333, 0.416822, 0.330333, 

....: -0.833644,-0.660666], [-1.0874,1.08734,-0.88595,0.885999],55.7898] 


Umjesto odeint, koristit čemo ode_solver, vidi Diferencijalne jednadžbe (alternative). 


sage: 
sage: 
sage: 
sage: 
sage: 
sage: 
sage: 


sage: 
sage: 


sage: 


S = ode_solver() 

S.function = system 

S.ode_solve{y_0 = figure8[ 0 ] , t_span=[0, figure8[2]], num_points = 1000) 
sol_figure8 = S.solution 

S.ode_solve{y_0=goggles[ 0 ], t_span=[0, goggles[2]], num_points=2000) 
sol_goggles = S.solution 

S.ode_solve{y_0=butterfly2[ 0 ], t_span=[0, butterfly2[ 2 ]], 

num_points=2000) 

sol_butterfly2 = S.solution 

S.ode_solve{y_0=yinyang2a[ 0] , t_span=[0, yinyang2a[ 2 ]], 

num_points=6000) 


sol_yinyang2a = S.solution 


Također, umjesto crtanja Matplotlibom koristit čemo Sageove grafičke elemente, kon kr etno primitivni 
grafički objekf line. 

sage: G = Graphics{) 

sage: xshift = 0.005 # for visibility 
sage: yshift = 0.01 # for visibility 

sage: G += line([{x[0],x[ 1] ) for t,x in sol_figure8]) 

sage: G += line([{x[2]+xshift,x[3]+yshift) for t,x in sol_figure8], 

. . . . : color='red') 
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sage 


sage 


G += line([{x[4]+2*xshift,X[5]+2*yshift) 
show(G, figsize = [6,6]) 


for t,x in sol_figureS], 

color='green') 



sage: 
sage: 
sage: 
sage: 
sage: 


G = Graphics {) 


G 

+ = 

line([{x[ 0 ] 

,x[l] ) 

for 

t , X 

in 

G 

+ = 

line([{x [2] 

,x[3] ) 

for 

t , X 

in 

G 

+ = 

line([{x [ 4 ] 

,x[5] ) 

for 

t , X 

in 


show(G, figsize = [6,6]) 


sol_goggles]) 
sol_goggles], 
sol_goggles], 


color=' red' ) 
color=' green' ) 



sage: G = Graphics{) 

sage: G += line([{x[0],x[ 1] ) for t,x in sol_butterfly2]) 

sage: G += line([{x[2],x [3] ) for t,x in sol_butterfly2], color='red') 


5.1. Mehanika 
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sage: G += line([{x[4],x[5]) for t,x in sol_butterfly2], color=' green' ) 
sage: show(G, figsize = [6,6]) 



sage: 
sage: 
sage: 
sage: 
sage: 


G = Graphics {) 

G += line([{x[ 0 ],x[ 1] ) for 
G += line([{x[2],x[3]) for 
G += line([{x[4],x[5] ) for 
show(G, figsize = [6,6]) 


t,x in sol_yinyang2a]) 
t,x in sol_yinyang2a], 
t,x in sol_yinyang2a], 


color=' red' ) 
color=' green' ) 



5.2 Elektrodinamika 
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Zadatak E-1 

Potencijal u točki r = (x, y, z) koji je posljedica statičke raspodjele N nabijenih čestica s 
nabojima gi, 92 , ■ ■ •, gAf i položajima fi = (xi, yi, zi), 7=2 = (x2, 2 / 2 , Z 2 ),..., tn = {xn, Un, zn) je 
dan formulom: 

N 

(t>ix,y,z) = 1^^^^ I • 

\r — ru\ 
k=i ' 

Definirajte funkciju plotpot (naboji, xgranice, ygranice) koja crta ekvipotenci- 
jalne konture u ravnini z = 0, gdje su naboji i njihovi položaji zadani kao lista oblika 
[(gi, xi, 2 / 1 , zi), (g 2 ,3^2, 2 / 2 , Z 2 ), ■ ■ .]. Upotrijebite tu funkciju da nacrtate ekvipotencijalne konture 
za naboje [ (-1, 1, 1, 1), ( + 1, -1, -1, 0.3), ( + 1, -1, 1, 0.5)]. 

Naputak: Koristite funkciju contour_plot {) da nacrtate cj){x, y, 0) 


5.2.1 Zračenje ubrzavanog naboja 

Snaga zračenja po prostornom kutu za točkasti naboj g brzine /? (u jedinicama brzine svjetlosti c) i 
ubrzanja a dana je Poyntingovim vektorom: 

dP |n X ((n —/3) X d)p 

dVL IbTT^eoc^ {I - fi ■ jdf 

gdje je n jedinični vektor u smjeru promatrača, a d jedinični bezdimenzionalni vektor ubrzanja. Cf. npr. 
Griffits,Eq. (11.72) 

sage: var('theta phi' ) 

(theta, phi) 


sage 


def unitn (theta, phi): 

"""Unit vector \hat{n} in Cartesian coordinate system.""" 

return vector((sin(theta)*cos(phi) , sin(theta)*sin(phi) , 

cos(theta))) 


Funkcija poynting() definira samo drugi bezdimenzionalni faktor u gornjem izrazu koji opisuje kutnu 
ovisnost 


sage 


def poynting (theta, phi, beta, ahat): 

"""Angular part of Poynting vector. 


nk = unitn(theta, phi) 

num = (nk.cross_product((nk-beta).cross_product(ahat))).norm()^2 
denom = (1-nk . dot_product (beta) ) "'S 
return (num/denom) * nk 


sage 


def poynting2d (theta, v, a): 

"""Projection of Poynting vector on phi=0, i.e x-z plane. 


P = poynting(theta, 0, v, a) 
return (p [ 0], p [ 2]) 


5.2. Elektrodinamika 
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Zračenje kad je akeeleraeija u smjeru brzine. Cf. Griffits, slika 11.14. 

sage: parametric_plot(poynting2d(theta, vector((0.9, 0, 0)), 

vector({l, 0, 0))), (theta, 0, 2*pi), axes_labels=['x', 'z']) 



Sinhrotronsko zračenje, kad je akeeleraeija okomita na brzinu. Cf. Griffits, slika 11.16. 

sage: parametric_plot(poynting2d(theta, vector((0.5, 0, 0)), 

vector((0, 0, 1))), (theta, 0, 2*pi), axes_labels=['x', 'z']) 



sage 


def plotField (v=vector(( 0. 5, 0, 0)), a=vector((l, 0, 0))): 

"""Plot the 3D shape of radiation field (cut along x-z plane) 


R = parametric_plot3d(poyntlng(theta, phi, v, a), (theta, 0, pi), 

(phi, 0, pl), plot_points=[100,100], frame=False, opacity=0.7) 
Av = arrow3d((0, 0, 0) , 3*v, color='green', wldth=5) 

Aa = arrow3d((0,0,0), 3*a, color='red', wldth=3) 

Tv = text3d("v", 3.6*v, slze=15) 

Ta = text3d("a", 3.3*a, slze=50) 
return R + Av + Aa +Tv + Ta 


sage: plotField () . show (inesh=True) 
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5.2.2 Ukupna snaga zračenja za sinhrotronsko zračenje: 


Za elektron u magnetskom polju B, okomitom na njegovu trenutnu brzinu, ubrzanje usljed Lorentzove 
sile je 

qcl3B 


a = 


rrip. 


Nakon uvrštavanja, rezultat je zgodno izraziti preko Thomsonovog udarnog presjeka 


(7t = 


2 /Uoe^ 


127rmgeoc^ 


= 6.652 m^ , 


i gustoće energije magnetskog polja 


UB = f- 

2/io 


pa se konačno dobije za ukupnu snagu zračenja: 

P = (jtcB^Ub / dn- -^L . 

47r ^ J (l-n-/3)5 

Iskoristiti ćemo ovaj primjer za demonstraciju računanja s mjernim jedinicama upotrebom SciPy modula 
constants i Pythonovog modula units. 


sage: from scipy import constants 

sage: c = constants.c * (units.length.metar / units.time.second) 
sage: q = units.charge.elementary_charge 

sage: sigmaT = constants.physical_constants [' Thomson cross section' 

. . . . : ] [0] * units . length.meter''2 

sage: muO = constants.mu_0 

sage: print "Thomsonov udarni presjek = " + str (sigmaT) 

Thomsonov udarni presjek = (6.652458734e-29)*meter^2 


5.2. Elektrodinamika 
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sage: UB{B) = B^2/2/constants.mu_0 * (units.energy.joule/units.length.meter^3) 

Relevantni integral po prostornom kutu: 

sage : def Kint (v=vector((0.5, 0, 0)), a=vector((0, 0, 1))): 

return integral_numerical(lambda theta: sin(theta) * 
integral_nuinerical (lambda phi : poynting (theta, phi, v, 

....: a).dot_product(unitn(theta, phi)), 0, 2*pi)[0], 0, pi)[0] 

sage: def pwr (beta, B): 

return (3/4/pi.n()) * sigmaT * c * beta''2 * UB(B) * Kint ( 

v=vector((beta,0,0))) 


sage: print "Ukupna snaga zračenja = " + str(pwr(0.5, 1)) 
Ukupna snaga zračenja = (7.05359483944320e-15)*joule/second 

Za kontrolu uspoređujemo s formulom 


P = 2aTf3^7^cUB 


sage: gamma(beta) = 1 / (l-beta"'2) 

sage: pwr2(beta, B) = 2 * sigmaT * beta^2 * gamma(beta) ^2 * c * UB(B) 
sage: pwr2(0.5, 1.) 

(7.0535948394432Oe-15)*joule/second 


5.3 Termodinamika i statistička fizika 

5.3.1 Brovvnovo gibanje 


Zadatak T-1 

Konstruirajte funkciju brown (n) koja daje listu [(xo, yo), {xi,yi),...] pozicija čestice koja se 
giba u ravnini slučajnim gibanjem koje je definirano rekurzijama Xi = Xi-i + r i yi = yi_i + s 
gdje su r i s slučajni brojevi između -1 i 1. Nacrtajte odgovarajuču putanju čestice. 

Naputak: 

• Inicijalizirajte listu pozicija tako da sadrži (0,0) kao prvi vektor pozicije i putem petlje 
dodajte u svakom koraku toj listi novi slučajni vektor. 
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Zadatak T-2 

Konstruirajte funkciju walklD {) tako da simulira tzv. jednodimenzionalnog nasumičnog šetača. 
Svaki šetačev korak treba biti iste, jedinične duljine, ulijevo ili udesno. Dakle, umjesto vektora 
(r, s) iz gornjeg Brovvnovog gibanja trebamo slučajan odabir koraka +1 ili -1. Uvjerite se u 
ispravnost tvrdnje da je prosječna kvadratna udaljenost nasumičnog jednodimenzionalnog šetača 
od ishodišta nakon N koraka N, sa standardnom devijacijom V^N, dok je prosječna apsolutna 
udaljenost 


5.4 Kvantna fizika 


5.4.1 Pravokutna potencijalna jama 

Isprogramirat ćemo funkciju koja računa kvantnomehaničke valne funkcije i pripadajuće energije veza¬ 
nih stanja čestice mase m u pravokutnoj potencijalnoj jami dubine V i širine L. Koristimo oznake s 
Wikipedijine stranice. 

sage: var ('a kmLVEABGH') 

(a, k, m, L, V, E, A, B, G, H) 


sage: import scipy.linalg as la 
sage: import numpy as np 
sage: import functools 


Definiramo valnu funkciju za potencijalnu jamu širine L. 


sage 

sage 

sage 

sage 

sage 

sage 


# components of piecewlse wave function 
psil(G, a, x) = G*exp{a*x) 

psi2(A, B, k, x) = A*sin(k*x) + B*cos(k*x) 
psi3(H, a, x) = H*exp{-a*x) 

# and a complete wave function 

def psi (x, coef s= ( 0 , 0 , 0 , 0 ) , wns={0, 0), L=l, norin=l, shift=0): 
""" coefs = (A, B, G, H) 
wns = (k, a) 

If !? !? 

if X < -L/2: 

psi = psil(coefs[2], wns[l], x) 
elif X > L/2: 

psi = psi3 (coefs [3], wns[l], x) 
else: 

psi = psi2(coefs[0], coefs[l], wns[0], x) 
return norm*psi + shift 


Valni brojevi u unutrašnjosti k i izvan jame a. 

sage: # hbar=l 

sage: wns = {k: sqrt (2*in*E) , a: sqrt (2*m* (V-E) ) } 


Rubni uvjeti na lijevom i desnom rubu jame su da tamo valna funkcija i njena prva derivacija budu 
neprekidne. 


5.4. Kvantna fizika 
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sage 


boundary_conditions = [ 


psil(G, a, - 

•L/2) == 

psi3 (H, a. 

L/2) == 

diff(psil(G, 

a, x) , 

diff(psi3(H, 

a, x) , 


] 


psi2 (A, B, k, -L/2), 
psi2 (A, B, k, L/2) , 

x) .subs(x=-L/2) == diff(psi2 (A, B, k, x) , 

x).subs(x=-L/2), 

x).subs(x= L/2) == diff(psi2 (A, B, k, x), 

x).subs(x= L/2) 


Ovi rubni uvjeti predstavljaju homogeni linearni sustav jednadžbi s nepoznanicama A, B, G i H. Energije 
vezanih stanja su određene zahtjevom da taj sustav ima rješenje. Konkretno, pretvorit ćemo rubne uvjete 
u matričnu jednadžbu oblika 


M \ 

M ^ =0 

Gr 

\H J 

a energije su onda dane kao rješenja obične jednadžbe det M = 0. Zbog trenutnog buga u Sageovoj ru¬ 
tini za determinantu matrica 4x4 i većih (trebao bi biti ispravljen u nadolazećoj verziji Sagea) definiramo 
svoju rutinu za izračun determinante: 

sage : def inydet (m) : 

"Evaluates determinant of 4x4 matrix m. (Sage's det{) has a bug.)" 

. . . . : d = 0 

. . . . : for c in range(4) : 

. . . . : ci = range (4) 

....: ci.pop(c) 

sub = m.inatrix_from_rows_and_columns { [1,2,3] , ci) 
d += (-l)^(c) * m[0,c] * det(sub) 

....: return d 


Elemente matrice M dobivamo metodom simboličkog izraza coef {), a jednadžbu detM = 0 rje¬ 
šavamo algoritmom u kojem krećemo od liste intervala [ (0, V) ], koja sadrži u početku samo jedan 
interval (0, V) i postupamo ovako: 

1. Maknemo prvi interval (a, b) iz liste (metoda pop {)) i tražimo u njemu nultočku. To može 
imati dva ishoda 

1. Našli smo nultočku (taj ishod je zagarantiran u prvom koraku jer jedno vezano stanje uvijek pos¬ 
toji). Zapišemo nađenu energiju, a listu intervala nadopunimo s dva nova intervala: (a, E) i 
(E, b). 

2. Nismo našli nultočku. Ne radimo ništa (interval je već maknut s popisa). 

2. Ponavljamo postupak iz točke 1. dok ne iscrpimo sve intervale. (Radi jednoznačnosti, svi intervali 
gore su u samom kodu zapravo za malu vrijednost eps udaljeni od navedenih rubova.) 

sage: def energies {system = {m: 1/2, L: 1, V: 25}, eps=le-3): 

bcs = [c.subs(wns).subs(system) for c in boundary_conditions] 

. . . . : # Preparing linear systein M*coefs = 0 

....: M = matrix([[(eq.lhs()-eq.rhs{)).coefficient(v) for v in (A,B,G,H)] 

. . . . : for eq in bcs]) 

. . . . : dt = mydet (M) 

. . . . : res = [ ] 

....: # Energies of bound States are ali Solutions of det(M) = 0 

. . . . : ints = [(eps, system[V]-eps)] 
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while ints: 


int = ints.pop(0) 


try: 

e = find_root(dt, int[0], int[l]) 
res.append(e) 

ints.extend{[(int[0], e-eps), (e+eps, int[l])]) 


except: 
pass 


res.sort() 
return res 


Za provjeru odredit ćemo energije i odgovarajuće vrijednosti bezdimenzionalne varijable v = A:L/2 za 
slučaj sa spomenute straniee na Wikipediji: 

sage: ens = energies(system = {m: 1/2, L: 1, V: 80}); ens 
[6.557920590133058, 25.767519828716722, 55.56613192980969] 

sage: # compare with viikipedia values = 1.28, 2.54, 3.73 

sage : [{k*L/2).subs(wns).subs{{m: 1/2, L: 1, V: 80}).subs(E=en).n() 

....: for en in ens] 

[1.28042186311124, 2.53808588451596, 3.72713468799458] 

Sada definiramo funkciju koja za datu energiju određuje koeficijente valnih funkcija A, B, G, H. Sustav 
jednadžbi rubnih uvjeta nije dovoljno odreden (detM = 0!). Dodatni uvjet koji bi potpuno odredio 
sustav je uvjet normalizacije valnih funkcija (integral gustoće vjerojatnosti po cijelom prostoru 

mora biti jedan), ali dodavanje tog uvjeta bi učinilo jednadžbe nelinearnima pa čemo raditi na način da 
privremeno fiksiramo normalizaciju valne funkcije uvjetom /2) = 1, što nam fiksira vrijednost 
koeficijenta H i onda za preostale koeficijente rješavamo nehomogenu linearnu matričnu jednadžbu 



sage: def coefs (en, system = [m: 1/2, L: 1, V: 25}): 


"""Returns (A, B, G, H) for given energy.""" 

# Fixing H by normalization psi(L/2)=1 

hval = numerical_approx(exp(a*L/2).subs(wns).subs(systein).subs( 


{E: 


en}) ) 


bcred = [c . subs (H=hval) . subs (wns) . subs (systein) . subs ({E : en}) 

for c in boundary_conditions] 
# Preparing linear systein Mred * coefs [:-l] - Fred = 0 
Mred = matrix([[(eq.Ihs()-eq.rhs()).coeff(v) for v in (A,B,G)] 

for eq in bcred]) 

Fred = vector([ (eq.Ihs()-eq.rhs()) .subs({A:0, B:0, G:0}) 

for eq in bcred]) 

cfs = la.lstsq(Mred, -Fred) [0] .tolist () 
cfs.append(hval) 
return tuple(cfs) 


Sada definiramo funkciju koja kombinira gornje funkcije, normalizira dobivene valne funkcije i vrača 
listu [( Ei, V’i), ( E 2 , 'ip 2 ), ■■■] gdje su funkcije jednog argumenta (x) dobivene parcijalnim 

izvrijednjavanjem (cf. currying ) na početku definirane funkcije psi za različite vrijednosti njenih 
opcionalnih argumenata, pomoću modula f unctools. (Normalizirane funkcije još množimo factorom 
scale radi čitlijvijeg crtanja.) 


5.4. Kvantna fizika 
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sage 


def sqwell ( systein = {m: 1/2, L: 1, V: 25}, scale=5.) : 
ens = energies(system) 
print "Energies = %s" % str (ens) 
res = [] 
for en in ens: 

cfs = coefs(en, system) 

#print cfs 

norm = l/sqrt(numerical_integral(lambda x: psi(x, cfs, 

(k.subs(wns) .subs(system) .subs (E=en), a.subs(wns) .subs ( 
systein) .subs (E=en) ), L=l)**2, (-oo, oo))[0]) 

#print en, 

fun = functools.partial (psi, coefs=cfs, 

wns=(k.subs(wns).subs({m: 1/2, L: 1, V: 25}).subs(E=en), 
a.subs(wns).subs(system).subs(E=en)), 

L=system[L], norin=norm*scale) 

#print fun(0) 
res.append((en, fun)) 
return res 


sage 


def square_well (system = {m: 1/2, L: 1, V: 25}, span=2, scale=5 . ): 
"""špan = horizontal špan of figure in widhts of well 

scale = scaling factor for wavefunctions for visibility 


colors = ['indigo', 'firebrick', ' darkolivegreen', 'plum'] 

sol = sqwell(system) 

P = line([[-span*system[L]/2, systein[V]], [-system[L]/2, 
system[V]], [-system[L]/2, 0], 

[system[L]/2, 0], [system[L]/2, system[V]], 

[ span*systein [L]/2, system[V]]], 
color='darkslategray', thickness=10, alpha=0.3) 
for (en, fun), k in zip(sol,range(len (sol))) : 

P += line([[-span*system[L]/2, en], [span*system[L]/2, en]], 

color=colors[k], thickness=0.4) 

P += plot(lambda x: fun (x, shift=en) , (-span*system[L]/2, 
+span*system[L]/2), color=colors[k]) 

P.show() 


sage: square_well(system = {m: 1/2, L: 1, V: 80}) 

Energies = [6.557920590133058, 25.767519828716722, 55.56613192980968] 
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5.5 Kozmologija 

Nobelova nagrada za fiziku 2011. godine dodijeljena je Perlmutteru, Schmidtu i Riessu “za otkriće 
ubrzanog širenja svemira putem opažanja dalekih supernova”. U ovom odjeljku ćemo prikazati taj 
rezultat. 

Udaljenost (luminosity distance) Dl {u megaparsecima Mpc) do supemove s prividnim sjajem m (mjeri 
se izravno) i apsolutnim sjajem M (poznat je za supernove tipa la) je dana formulom p = m — M = 
5 logiQ Dl + 25. Na internet adresi http://supernova.lbl.gov/Union/ nalaze se podaci o nekoliko stotina 
supemovih organiziranih u tablici sa stupcima: 

1 . ime supemove 

2 . z - crveni pomak 

3. p - atenuacija sjaja supemove 

4. Ap - eksperimentalna neodređenost (greška) od p 

Učitavamo ih koristeći NumPy funkciju loadtxt {) koja numeričke podatke organizirane u stupce 
obične datoteke pretvara u NumPy listu. 


sage : 

import matplotlib.pyplot as pit 


sage : 

import numpy as np 



sage : 

import scipy 



sage : 

import scipy.stats 



sage : 

import urllib 



sage : 

sns = urllib.urlopen( 




'http://supernova.Ibl 

.gov/Union/figures/SCPUnion2_mu_vs 

sage: 

dat = np.loadtxt(sns, 

usecols= (1,2,3) 

; dat.shape 

(557, 

3) 



sage : 

dat [ : 2 , : ] 



array ( [[ 2.84880000e-02, 

3.533555110+01, 

2.261439260-01], 


[ 5.004300000-02, 

3.667544160+01, 

1.671142520-01]] 
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Pretvaramo sada mjerene podatke iz ^ i A/r u podatke za Dl i ^Dl pazeći na ispravnu propagaciju 
greške. 

sage: var ('mu DL') 

(mu, DL) 

sage: DL(mu) = DL.subs(solve(mu == 5*log(DL, base=10) + 25, DL)[0]); DL 
mu I—> 1/100000*0''(l/5*mu*log (10) ) 


sage: DLdat 


(557, 3) 


np.array ( [ (z, DL(valmu) .n(), diff(DL) (valmu) .n() *errmu) 

for z, valmu,errmu in dat]); DLdat.shape 


Definiramo funkcije za prilagodbu, poznate iz odjeljka o prilagodbi funkcija, a dodajemo i dio za ocjenu 
dobrote prilagodba, vidi odjeljak o testiranju hipoteze. 


sage 


sage 


sage 


def dist (p, fun, data) : 

return np.array([(y - fun(p, x))/err for (x,y,err) in data]) 

def chisq(p, fun, data): 

return sum ( dist (p, fun, data)''2 ) 

def leastsgfit (fen, init, data): 

"""Fit funetion fen, starting from init initial values to data.""" 
p_final, eov_x, infodiet, msg, ier = seipy.optimize.minpaek.leastsg( 

dist, init, (fen, data), full_output=True) 
parerrs = sqrt(seipy.diagonal(eov_x)) 
for k in range(len(init)) : 

print "p[%d] = %.3f H— %.3f" % (k, p_final[k], parerrs[k]) 
ehi2 = ehisq(p_final, fen, data) 
df = len(data)-len(init) 

print " ehi''2/d. o. f = %.lf/%i (p-value = %.4f)" % (ehi2, df, 
seipy.stats.ehisgprob(ehi2, df) ) 
return list(p_final) 


5.5.1 Određivanje Hq iz podataka o biiskim supernovama 

Za male crvene pomake, ovisnost udaljenosti i crvenog pomaka dana je Hubbleovim zakonom 


HoDl = cz 


gdje je c brzina svjetlosti, a Hq Hubbleova konstanta. Određujemo Hubbleovu konstantu (u jedinicama 
km/s/Mpc) prilagodbom Hubbleovog zakona na podskup mjerenja za koja je t < 0.12. 

sage: lowz = np.array([pt for pt in DLdat if pt[0]<0 . 12]); lowz.shape 
(174, 3) 

sage: e=3.e5 

sage: def linfun (p, z) : 

....: return c*z/p[0] 

sage: HO, = leastsgfit(linfun, (60,), lowz) 
p[0] = 68.012 +- 0.385 

chi^2/d.o.f = 179.1/173 (p-value = 0.3585) 

sage: print "\nHubbleova constanta HO = %.lf km/s/Mpc" % HO 
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Hubbleova constanta HO = 68.0 km/s/Mpc 


sage: fig, ax = pit.subplots{figsize=[5,3]) 

sage : ax . errorbar {lowz [ :, 0 ] , lowz yerr=lowz [ :, 2 ] , marker=' .' , 

. . . . : linestyle='None') 

sage: zvals = np.linspace (0 . 01, 0.12) 

sage: ax.plot(zvals, [linfun([HO], z) for z in zvals], 'r-') 
sage: ax.set_xlabel( '$z$' ) 

sage: ax.set_ylabel( '$D_L \; {\\rm [Mpc]}$') 

sage: fig.tight_layout() 
sage: fig.savefig{ ' fig' ) 



5.5.2 Određivanje ubrzanja ekspanzije svemira 

Za veće ervene pomake, Hubbleov zakon se modifieira i u prvoj aproksimaeiji se može pisati kao 

HqDl = c 2;(1 + (1 - go)|) 

gdje je go parametar deeeleraeije. Definiran je tako da pozitivni go odgovara svemiru čija se brzina 
širenja usporava (kako se očekivalo prije ovih mjerenja). Određujemo taj parametar prilagodbom ove 
formule na eijeli skup mjerenja. 

sage: def qfun (p, z) : 

....: return c*z/p[0] * (1 + 0.5*{1 - p[l])*z) 

sage: gpars = leastsgfit(qfun, (50, 0.1), DLdat) 
p[0] = 69.184 +- 0.368 
p[l] = -0.153 +- 0.039 

chi^2/d.o.f = 553.7/555 (p-value = 0.5081) 


sage: fig, ax = pit . subplots(figsize=[5, 3] ) 

sage: ax . errorbar(DLdat[:, 0 ], DLdat[:, 1] /le3, yerr=DLdat[:, 2 ]/le3, 

. . . . : marker=' .', linestyle='None') 

sage: zvals = np . linspace ( 0 . 01 , 1.8) 


5.5. Kozmologija 
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sage: ax.plot(zvals, qfun(qpars, zvals)/le3, 'r' ) 
sage: ax.set_xlabel( 'Sz$' ) 

sage: ax.set_ylabel( '$D_L \; {\\rm [Gpc]}S') 

sage: fig.tight_layout() 
sage: fig.savefig{ ' fig' ) 





Iako iz mjernih podataka nije očito da postoji odstupanje od linearnog Hubbleovog zakona, podataka 
ima dovoljno da statistička analiza nedvosmisleno pokazuje da je parametar deceleracije negativan i 
različit od nule: 


qo = -0.15 ±0.04 


5.5.3 Određivanje gustoće materije i tamne energije 

Gornja formula ne vrijedi za z « 1, dakle ovakav račun nije sasvim korektan i može se rabiti samo 
u pedagoške svrhe. Ispravan pristup zahtijeva integraciju propagacije zrake svjetlosti od supernove do 
promattača kroz svemir u širenju. Ta propagacija je određena sastavom svemira (obična i tamna materija 
različito djeluju od tzv. tamne energije koja uzrokuje ubrzanje širenja svemira). Hubbleov zakon postaje 

HoDl = c(l ± z)Z{z, Qm, ^a) 


Z(^z^ j ^a) 


da 


'i/{i+z) a^^X{a,VlM,^A) 


X{a, ^Mi ^ a ) 


a 




rijvf ± — 1 

gdje je fŽM udio materije (obične i tamne), a Oa udio tamne energije u ukupnoj energiji svemira.Još 
općenitije formule, koje ne rade pretpostavku Qm ± ^^A = 1, mogu se naći na stranicama Neda Wrighta. 

Provjeravamo egzaktne formule razvojem u red po z: 
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sage: var('z a Om Ov Z X q0' ); assume(z>0, a>0, Om>0, Ov>0); 

(z, a. Om, Ov, Z, X, q0) 

sage: X{a, Om, Ov) = Om/a + Ov*a**2 + (l-Om-Ov) 

sage: DLHO = taylor((( 1 + z) ^2 * (Z/(1 + z)) * (1 +(1-Om-Ov)*Z^2/6)) .subs ( 
....: Z==taylor(integral(1/(a*sqrt(X(a. Om, Ov))), 

a, l/(l+z), 1), z, 0, 2)), z, 0, 2) 

sage: solve(DLHO == z*(l + (1 - q0)*z/2), q0)[0] 

q0 == 1/2*Om - Ov 


Ova korespondencija s parametrom deceleracije se slaže s (Peacock 3.34). 

sage: def J (x): 

if x<0: 

return sin(sqrt(-x))/sqrt(-x) 
elif x==0: 

return 1 
else: 

return sinh(sqrt(x))/sqrt(x) 


sage: def Xfun(a, Om, Ov): 

....: return Om/a + Ov*a**2 + (1-Om-Ov) 


sage: def Zfun(z, Om, Ov): 

return numerical_integral(lambda a: 1/(a*sqrt(X(a. Om, Ov))), 

. . . . : 1/(1 + z) , 1) [0] 

sage: def DLfun (p, z) : 

....: """DL(z) with pars p=[H0, Omega_m, Omega_Lambda] 

. . . . : HO = p[0] 

. . . . : Om = p [ 1] 

. . . . : Ov = p [ 2] 

. . . . : Otot = Om + Ov 

. . . . : Zval = Zfun (z. Om, Ov) 

return c/HO * (1+z) * Zval * J((1-Otot)*Zval**2) 


sage: def DLfunflat (p, z) : 

....: """DL(z) with pars p=[H0, Omega_Lambda]. Fiat universe. 

. . . . : HO = p[0] 

. . . . : Ov = p[1] 

. . . . : Om = 1-Ov 

. . . . : Zval = Zfun(z. Om, Ov) 

....: return c/HO * (1+z) * Zval 


Prilagodbom parametara Hq, Q.m i određujemo udjele pojedinih komponenti svemira 

sage: fullpars = leastsgfit(DLfun, (70, 0.3, 0.), DLdat) # long time 
p[0] = 70.717 +- 0.466 
p[l] = 0.368 +- 0.075 
p[2] = 0.834 +- 0.122 

chi^2/d.o.f = 527.7/554 (p-value = 0.7833) 

Ukoliko pak fiksiramo ravnu geometriju ( Qm + = 1) dobivamo uobičajenu jednu trećinu materije 

i dvije trećine tamne energije: 
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sage: flatpars = leastsqfit(DLfunflat, (70, 0.3), DLdat) # long time 

sage: print "Om = 1-Ov = %8.3f" % {1-flatpars [1] , ) 
p[0] = 70.412 +- 0.363 
p[l] = 0.708 +- 0.021 

chi"2/d.o.f = 528.7/555 (p-value = 0.7827) 

Om = 1-Ov = 0.292 

Dakle, 


VLm — 0.29 
Da = 0.71 

Crtamo ovaj rezultat (crvena linija) zajedno s nekim drugim karakterističnim scenarijima. 

sage: fig, ax = pit.subplots{figsize=[8, 6] ) 

sage: ax.errorbar(DLdat[:,0], DLdat[:,1], yerr=DLdat[:,2], marker=' .' , 

. . . . : linestyle='None') 

sage: zvals = np.linspace(0 . 01, 1.8) 

sage: ax.plot(zvals, [DLfunflat(flatpars, z) for z in zvals], 'r-', 

....: label='$\\Omega_m=0.29, \\Omega_\\Lambda=0.71$') 

sage: ax.plot(zvals, [DLfun([H0, 0.5, 0.5], z) for z in zvals], 'b:', 
....: label='$\\Omega_m=0.5, \\Omega_\\Lambda=0.5$') 

sage: ax.plot(zvals, [DLfun([H0, 1, 0], z) for z in zvals],' k-.' , 

....: label='$\\Omega_m=l, \\Omega_\\Lambda=0$') 

sage: ax.plot(zvals, [DLfun([H0, 0, 1], z) for z in zvals], 'g— .', 

....: label='$\\Omega_m=0, \\Omega_\\Lambda=l$') 

sage: ax.plot(zvals, [DLfun([H0, 0.3, 0.], z) for z in zvals], 'k— ', 
....: label='$\\Omega_m=0.3, \\Omega_\\Lambda=0$') 

sage: ax.set_ylim(0, 1.5e4) 
sage: ax.set_xlabel( '$z$' ) 

sage: ax.set_ylabel( '$D_L \; {\\rm [Mpc]}$') 

sage: ax.legend(loc=' upper left' ).draw_frame(0) 
sage: fig.tight_layout() 
sage: fig.savefig( 'fig' ) 
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POGLAVLJE 6 


Još matematike 


6.1 Kompleksna analiza 

Matematičke funkcije se redovito ispravno ponašaju i za kompleksne vrijednosti argumenta 

sage: log(-l+le-8*I) 

3.14159264358979*1 
sage: log(-l-le-8*I) 

-3.14159264358979*1 

Vidimo da logaritamska funkcija uredno radi s kompleksnim argumentom i daje nam glavnu granu 
multifunkcije s argumentom —vr < arg{log{z)) < vr, s diskontinuitetom od 27r prilikom prelaska 
negativne realne osi. 

sage : var( ' x, y' ) 

(x, y) 

sage: plot3d(imag(log (x+y* 1 j )), (x,-2,2), {y,-2,2)) 
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taylor () može dati i Laurentov razvoj u kompleksnoj ravnini. Npr, vrijedi 



što do nultog reda dobivamo ovako: 

sage : var( ' z' ) 
z 

sage: taylor(l/(z^2 + l), z. I, 0) 
-1/(2*z - 2*1) +1/4 


Za izračun reziduuma ne postoji posebna funkcija, ali i njega možemo dobiti pomoću funkcije 
taylor () zahvaljujući činjenici daje reziduum jednak minus prvom koeficijentu (onom uz l/{z — zo)) 
u Laurentovom redu. Npr. tražimo reziduum funkcije f{z) = exp{z)(z^ u z = 0: 

sage: Irnt = taylor(e^z/z^5, z, 0, -1); Irnt 
1/24/z + l/6/z^2 + l/2/z-'3 + l/z''4 + l/z''5 


Možemo koristiti i formulu s limesom u polu: 

sage: [limit (diff {z''k*e^z/z''5, z, k-1), z=0)/factorial (k-1) 
[+Infinity, -Infinity, +Infinity, -Infinity, 1/24] 


for k in 

range(1,6)] 


Obje metode ispravno daju 1/24 za traženu vrijednost reziduuma. 
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Zadatak 1 

Definirajte funkciju plotpoint (z) koja crta kompleksnu ravninu sa točkom koja odgovara 
kompleksnom broju z. Dakle x=Re(z), y=Im(z) 


Zadatak 2 

Definirajte funkciju argand (lista) koja crta kompleksnu ravninu s točkama zadanim u listi 
kompleksnih brojeva lista = [zl, z2, ...]. Upotrijebite je da nacrtate na jednom dija¬ 

gramu svih devet devetih korijena od 1. 


Zadatak 3 

Definirajte funkciju plot f un (f un, min, max ) koja crta realni i imaginarni dio funkcije f un 
za interval apscise (min, max). 


6.2 Vektorska polja 

Neka se čestica u ravnini nalazi u potencijalu U{x,y) = —Ax^y‘^. Odredimo odgovarajuču silu 
F = —Vf7 To je moguče pomoču metode simboličkih izraza gradient () koja vrača odgovarajuče 
vektorsko polje sile: 

sage : var( ' x y z' ) 

(x, y, z) 

sage: U = -4*x^2*y"'3 
sage: F = -U.gradient(); F 
(8*x*y^3, 12*x''2*y''2) 

Dakle, polje sile je F = i+12x^y^j i možemo ga skicirati funkcijom plot_vector_f ield {) : 

sage: plot_vector_field(F, (x,-2,2), (y,-2,2)) 



Ukoliko potencijal ovisi i o nekim parametrima moramo eksplicitno deklarirati varijable koordinatnog 
sustava. Isto trebamo i kad potencijal ne ovisi o nekoj koordinati. Npr. ukoliko je gornji potencijal 
zapravo U (x, y, z), dakle u 3D prostoru samo što ne ovisi o z, onda imamo: 


6.2. Vektorska polja 
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sage: F3 = -U.gradient([x,y,z]); F3 
(8*x*y^3, 12*x^2*y''2, 0) 

Sage nema specijalizirane funkcije za divergenciju i rotaciju, ali nije ih teško napisati. Jedan primjer 
nalazi se u datoteci dostupnoj ovako: 

sage : load ( ' http : / /www. phy . hr/~kk:uiner/sp/divcurl. py' ) 

- Loaded functions: div, curi 


sage: div(F3, [x,y,z]) 
24*x^2*y + 8*y^3 
sage: curl(F3, [x,y,z]) 
( 0 , 0 , 0 ) 


Provjerimo daje divergencija rotacije proizvoljne funkcije nula: 


sage: div(curi([x^2 + 2*y, x^3 + 

0 

Crtanje 3D vektorskog polja: 

sage: plot_vector_field3d(F3, (x, 

2.0 



0.0 


in(z), y*z + 2], [x,y,z]), [x,y,z]) 

2,2), {y,-2,2), (z,-2,2)) 

/ 

i 

i 



2.6-2.0 


6.3 Testiranje hipoteze 


često je potrebno kvantificirati dobrotu neke prilagodbe te odrediti da li predloženi teorijski model dobro 
opisuje mjerenja. Uzmimo slijedeći niz mjerenja s greškama (xj, ?/j, Oi = At/*): 

sage: errdata = [[0.55, 0.74, 0.20], [0.9, 1.07, 0.12], 

....: [1.2, 1.28, 0.12], [1.5, 1.02, 0.12], [1.9, 1.63, 0.20]] 

sage: xs = [x for x,y,err in errdata] 
sage: ys = [y for x,y,err in errdata] 
sage: errs = [err for x,y,err in errdata] 
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Za crtanje točaka s greškama (errorbars), koristimo Matplotlib funkciju errorbar () gdje greške idu 
u opcionalni argument yerr: 

sage: import matplotlib.pyplot as pit 
sage: import n\ampy as np 

sage: fig, ax = pit.subplots{figsize=[5,3]) 

sage: ax.errorbar(xs, ys, Yerr=errs, linestyle=' None' ) 
sage: fig.savefig{ 'fig' ) 



Promotrimo tri moguća teorijska opisa ovih mjerenja: potenciju logaritam log(ax) i obični pravac 
a + bx i odredimo parametre prilagodbom pomoću funkcije find^t(). (Zasad zanemarimo različite 


greške pojedinih mjerenja.) 




sage: 

(a, b) 

var( 'a, b' ) 




sage: 

powmodel (x) = x''a 




sage: 

logmodel(x) = log(a*x) 



sage: 

linmodel(x) = a*x+b 




sage: 

powfit=find_fit([[x. 

Y] 

for 

x,y,err in errdata], powmodel, 
solution_dict=True); print powfit 

{a : 0 . 

6166290268412898} 



sage: 

logfit=find_fit([[x. 

Y] 

for 

x,y,err in errdata], logmodel, 
solution_dict=True) ; print logfit 

{a: 2. 

836898540272243} 



sage: 

linfit=find_fit([[x. 

Y] 

for 

x,y,err in errdata], linmodel. 





solution dict=True); linfit 


{b: 0.49690475972372317, a: 0.5380952396683008} 


sage: 
sage: 
sage: 
sage: 


sage: 


xvals = np.linspace(0 .4 , 2) 

ax.plot(xvals, [powmodel(x).subs{powfit).n( ) for x in xvals] 
ax.plot(xvals, [logmodel(x).subs{logfit).n() for x in xvals] 
ax.plot(xvals, [linmodel(x).subs{linfit).n() for x in xvals] 

color='black', linestyle=' 


fig.savefig {' fig' ) 


,'r') 


'b-.' ) 
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Teško je “od oka” procijeniti koji je opis najbolji, a još teže koji opisi su prihvatljivi a koji nisu. Stan¬ 
dardna mjera dobrote fita je 


X 


2 


{Vi - 

2 -^ ^2 


gdje su cTj greške mjerenja. (Usput, kad je riječ o mjerenju frekvencija onda je greška dana kao ai = 
Dobre su prilagodbe s ~ df, gdje je df broj stupnjeva slobode (broj mjerenja minus broj 
slobodnih parametara funkcije koju prilagodavamo). Za točnije kvantificiranje dobrote prilagodbe gleda 
se integral statističke raspodjele s df stupnjeva slobode od izračunatog do beskonačnosti koji se 
često zove p-vrijednost i koji ne bi trebao biti manji od 0.1 ili barem 0.01 ili hipotezu treba odbaciti. 
Taj integral je implementiran kao funkcija scipy . stats . chisgprob (chisq, df ). (Inače, za 
razliku od ovog primjera gdje tražimo statističku podršku hipotezi da funkcije opisuju mjerenja, često 
smo u obrnutoj situaciji gdje testiramo tzv. nul-hipotezu. Tu se npr. pitamo koja je vjerojatnost nekih 
mjerenja ukoliko neka čestica ne postoji ili ukoliko neki neki lijek ne djeljuje. Tada je mala p-vrijednost 
“poželjna” jer predstavlja statističku podršku postojanju te čestice ili dobrom djelovanju lijeka.) 

Definirajmo tri funkcije koje prilagodavamo ... 


sage 


sage 


sage 


def powfun (p, x): 
return x''p[0] 

def logfun (p, x): 

return log(p[0]*x) 

def 1 infun (p, x): 

return p[l]*x+p[0] 


... i odgovarajuću funkciju udaljenosti tj. odstupanja (čije kvadrate će minimizirati leastsg ()). Tu 
trebamo staviti i težinski faktor l/ai za svaku točku kako bi točke s manjom greškom više utjecale na 
prilagodbu. Ovdje ćemo dodati još jedan argument putem kojeg ćemo prosljeđivati jednu od gornje tri 
funkcije: 

sage: def dist (p, fun, data): 

return np.array([{y - fun (p, x))/err for (x,y,err) in data]) 
Definirajmo i odgovarajući x^- 
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sage: def chisq(p, fun, data): 

....: return sum( dist(p, fun, data )^2 ) 

Nakon minimizacije ćemo testirati i zastavicu ier tako da ukoliko least sq () nije uspješna dobijemo 
odgovarajuću poruku o grešci: 

sage: import scipy 
sage: import scipy.stats 
sage: fitfun = powfun 

sage: ppow_final, cov_x, infodict, msg, ier = scipy.optimize.minpack.leastsq( 

dist, [1.], (fitfun, errdata), full_output=True) 

sage: if ier not in (1,2,3,4) : 

. . . . : print "-> No fit! <-" 

....: print msg 

....: else: 

print "a =%8.4f H— %.4f" % (ppow_final, sgrt(cov_x[0,0])) 

....: chi2 = chisq([ppow_final], fitfun, errdata) 

....: df = len(errdata)-1 

print "chi-'2 = %.3f" % chi2 
....: print "df = %i" % df 

print "p-vrijednost = %.4f" % scipy.stats.chisgprob(chi2, df) 

a = 0.5055 +- 0.1486 

chi"2 = 7.880 
df = 4 

p-vrijednost = 0.0961 


sage: fitfun = logfun 

sage: plog_final, cov_x, infodict, msg, ier = scipy.optimize.minpack.leastsg( 
....: dist, [1.], (fitfun, errdata), full_output=True) 

sage: if ier not in (1,2,3,4) : 

. . . . : print "-> No fit! <-" 

....: print msg 

....: else: 

....: print "a =%8.4f +- %.4f" % (plog_final, sgrt(cov_x[0,0])) 

....: chi2 = chisg([plog_final], fitfun, errdata) 

....: df = len(errdata)-1 

....: print "chi-'2 = %.3f" % chi2 

....: print "df = %i" % df 

....: print "p-vrijednost = %.4f" % scipy.stats.chisgprob(chi2, df) 

a = 2.7219 +- 0.1693 

chi^2 = 15.973 
df = 4 

p-vrijednost = 0.0031 


sage: fitfun = linfun 

sage: plin_final, cov_x, infodict, msg, ier = scipy.optimize.minpack.leastsg( 
....: dist, [1., 1.], (fitfun, errdata), full_output=True) 

sage: if ier not in (1,2,3,4) : 

. . . . : print "-> No fit! <-" 

....: print msg 

....: else: 

....: parerrs = sgrt(scipy.diagonal(cov_x)) 

....: print "p[0] =%8.3f +- %.4f" % (plin_final[0], parerrs[0]) 

....: print "p[l] =%8.3f +- %.4f" % (plin_final[1], parerrs[l]) 

....: chi2 = chisg(plin_final, fitfun, errdata) 

....: df = len(errdata)-2 # dva parametra 
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print "chi^2 = %.3f" % chi2 
....: print "df = %i" % df 

print "p-vrijednost = %.4f" % scipy.stats.chisgprob(chi2, df) 

p[0] = 0.656 +- 0.2121 

p[l] = 0.398 +- 0.1683 

chi"2 = 7.116 
df = 3 

p-vrijednost = 0.0683 


sage: fig, ax = pit.subplots{figsize=[5,3]) 

sage: ax.errorbar(xs, ys, yerr=errs, linestyle=' None' ) 

sage: xvals = np.linspace ( 0.4 , 2) 

sage: ax.plot(xvals, [powfun([ppow_final], x) for x in xvals],'r') 
sage: ax.plot(xvals, [logfun { [plog_final ] , x) for x in xvals ] , ' b-. ' ) 
sage: ax.plot(xvals, [linfun{plin_final, x) for x in xvals], color=' black' , 
. . . . : linestyle=' :') 

sage: fig.savefig{ ' fig' ) 



Ova analiza sugerira da možemo odbaciti logaritamski model, dok su druga dva modela, u nedostaku 
boljih, prihvatljiva. 

— 

Zadatak 4 

Prilagodite funkciju f{x) = ax donjim podacima. Odredite parametar a i njegovu grešku te ocije¬ 
nite dobrotu fita. 

sage: xs2, ys2, errs2 = [range(2,25,2) , [5.3, 14.4,20.7,30.1,35.0, 41.3, 52.7, 

....: 55.7,63.0,72.1,80.5,87.9], [1.5 for k in range(2,25,2) 

sage: errdata2 = zip(xs2, ys2, errs2) 
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POGLAVLJE 7 


Još programiranja 


7.1 Mogućnosti ispisa rezultata: 


Sage notebook sučelje može ispisati rezultat računa u različitim formatima: 

1. Defaultni ispis jedini omogućuje izravni copy/paste u input ćelije. 

sage: var('x n theta' ) 

(x, n, theta) 

sage: si =suin (sin {theta+x)/log (n+x) , x, 1 , oo) 

sage: si 

sum(sin(theta + x)/log(n + x), x, 1, +lnfinity) 

sage: s2 = matrix(ZZ, [[ 1 , 2, 3], [4, 5, 6]]) 

sage: s2 


[1 2 3] 
[4 5 6] 


2. Ljepši ispis koji koristi interni LaTeX i tzv. jsMath dobiva se funkcijom show (), ili aktiviranjem 
“Typeset” gumba na vrhu notebook sučelja. 

sage: show(sl) 
sage: show(s2) 




3. Ispis odgovarajućeg koda za copy/paste u LaTeX dokument. 

sage: latex(s2) 

\left ( \begin {array}{rrr} 

1 & 2 & 3 \\ 

4 & 5 & 6 
\end{array} \right) 

4. Prikaz LaTeX->PDF inačice u posebnom prozoru (Moguće je da ovo radi samo ukoliko su Sage 
klijent i server na istoj mašini ili uz pažljivo podešene X-windows autorizacije.) 


83 







Sage računalno okruženje za fizičare, Distribucija 1.1 


sage: view(s2, viewer=' pdf' ) 


7.2 Još o funkcijama 


Funkciju definiranu po dijelovima možemo dobiti na slijedeći način: 


sage 


sage: 


def fpw (x): 
if x<-l: 

return 
elif x>l : 

return 
else: 

return 

plot(fpw, -2 , 


-X 

X 

x -'2 

2 ) 



Python funkcije ne možemo općenito npr. simbolički integrirati. Posebno je opasno to što, ukoliko 
pokušamo, možemo dobiti pogrešan odgovor: 

sage: integral(fpw{x), x, 1, 2) 

7/3 

No, uvijek možemo pribjeći numeričkoj integraciji koja u ovom slučaju daje točan rezultat 3/2: 

sage: numerical_integral(fpw, 1, 2) [0] 

1.5 

Pri gornjem pokušaju simboličke integracije, prvo je izvrijednjen sam integrand fpw (x) . Kako u tom 
trenutku x nema nikakvu vrijednost, on ne zadovoljava ni jedan od dva uvjeta (x> 1 ili x< 1) pa funkcija 
vrača simbolički izraz x''2 (is else bloka), što je pogrešno za integraciju u granicama od 1 do 2. 

sage: fpw(x) 
x^2 

Iz istog razloga je prilikom crtanja gore bilo potrebno kao prvi argument staviti samu funkciju fpw, a 
ne izraz fpw (x) koji bi dao krivi crtež: 

sage: plot(fpw(x), x, -2, 2) 
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Naravno, ukoliko znamo da python funkcija uvijek vraća korektan simbolički izraz, smijemo je integri¬ 
rati i diferencirati: 

sage: def gun(f, x) : 

. . . . : return f{f(x)) 

sage: diff(gun (sin, x), x) 
cos(x)*cos(sin(x)) 


7.3 Funkcionalno programiranje 

Isti matematički problem je često moguće riješiti na različite načine, putem algoritama iza kojih stoje 
sasvim različiti načini razmišljanja. Razložimo to na primjeru funkcije faktorijel. 

sage: factorial(7) 

5040 

Klasični način programiranja, upotrebljavan od dana prvih računala, je tzv. proceduralno (ili impera¬ 
tivno) programiranje kod kojeg na program gledamo kao na niz naredbi koje obično postupno mijenjaju 
vrijednosti nekih varijabli pohranjenih u memoriji računala: 

sage: def factorialProc (n): 

. . . . : fac = 1 

. . . . : i = 1 

....: while i<n: 

. . . . : i = i + 1 

....: fac = fac * i 

....: return fac 

sage: factorialProc (7) 

5040 

Ovdje je očito riječ o standardnom algoritmu kakvog bismo isprogramirali u bilo kojem proceduralnom 
jeziku poput Fortrana ili C-a: Imamo petlju koja se prolazi n puta i svaki puta množimo rezultat prošlog 
prolaza sa za 1 večim brojem. 

Međutim, postoje i drugačiji pristupi. Tako je kod funkcionalnog programiranja naglasak na izvrijed- 
njavanju funkcija tj. primjeni funkcija na izraze. Faktorijel prirodno zamišljamo kao funkciju koja daje 
“umnožak svih brojeva od 1 do zadanog broja”. 

Kao prvo, treba nam modul operator koji implementira funkcije koje odgovaraju standardnim mate¬ 
matičkim operacijama *, +, -,... 
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sage: import operator 

sage: (operator.mul(2,3), operator.add{ 2 , 3 )) 

(6, 5) 

Zatim koristimo funkciju reduce {), koja kumulativno primjenjuje prvi argument (koji mora biti funk¬ 
cija) na elemente drugog argumenta: 

sage: f = function { 'f' ) 
sage: reduce (f, range(l,5)) 
f(f(f(l, 2), 3), 4) 


sage: reduce (operator.mul, range(l,8)) 

5040 

Primijetite da nam ovdje nisu trebale pomoćne varijable (poput i i f ac iz f actorialProc), ali nam 
je trebalo više memorije jer smo trebali pohraniti čitavu listu [ 1, 2 , . . . ] koju daje range (). 

Funkcionalno programiranje je stil programiranja koji stavlja naglasak na izvrijednjavanje izraza, a ne na 
sukcesivno izvršavanje komandi. Kod čistih funkcionalnih programa obično nema pomočnih varijabli 
i nepotrebnih popratnih efekata izvrijednjavanja funkcija. U tom smislu funkcionalno programiranje je 
pomalo slično radu s tabličnim kalkulatorom (npr. MS Excel): redosljed izračunavanja čelija je nebitan 
tj. očekujemo automatsku konzistenciju. Isto, vrijednosti ćelija su dane izrazima, a ne slijedovima 
komandi. (Misli se na Excel, a ne Sage ćelije.) 

(Više informacija o funkcionalnom programiranju nudi Eunctional Programming PAQ ili WWW stranice 
Haskell funkcionalnog jezika.) 

Takav stil programiranja često rabi nekoliko specijalnih funkcija (mogli bismo ih zvati “meta-funkcije”) 
čija je uloga upravljanje primjenivanjem drugih funkcija koje dolaze kao argumenti ovih meta-funkcija. 
Možda najvažnija takva funkcija je map () koja distribuira (mapira) funkciju po elementima neke liste: 

sage: f = function (' f' ) 
sage: map(f, range (4)) 

[f(0), f(l), f(2), f(3)] 

Ukoliko funkcija prima više argumenata, može se mapirati na više listi: 

sage: map(f, range (4), range (3,7)) 

[f(0, 3), f(l, 4), f(2, 5), f(3, 6)] 

Recimo da sad želimo ispisati tablicu s nizom prirodnih brojeva i njihovih faktorijela. Možemo prvo 
postupiti tako da prvo definiramo pomoćnu funkciju koja generira jedan red tablice i onda je pomoću 
map () primijenimo na niz brojeva: 

sage: def auxfun(k:) : 

. . . . : return (k, factorial (k)) 


sage: map(auxfun, range (7)) 

[(0, 1), (1, 1), (2, 2), (3, 6), (4, 24), (5, 120), (6, 720)] 

Međutim, baš kao i pomoćne varijable, tako ni pomoćne funkcije nisu u duhu funkcionalnog programi¬ 
ranja. Zato postoje tzv. lambda-funkcije koje su bezimene i definiraju se i upofrebljavaju na slijedeći 
način: 
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sage: map(lambda k: (k, factorial(k)) , range(7)) 

[(0, 1), (1, 1), (2, 2), (3, 6), (4, 24), (5, 120), (6, 720)] 


sage: matrix(_) # čitljiviji ispis bez uptrebe formatirajućih stringova 
[ 0 1 ] 

[ 1 1 ] 

[ 2 2 ] 

[ 3 6] 

[ 4 24] 

[ 5 120] 

[ 6 720] 


Treba uočiti kako je često upotrebu funkcije map () i lambda-funkcija moguče izbječi korištenjem obu¬ 
hvaćanja liste. Npr, gornji primjer se elegantnije realizira ovako: 

sage: [ (k, factorial(k)) for k in range{7)] 

[(0, 1), (1, 1), (2, 2), (3, 6), (4, 24), (5, 120), (6, 720)] 


Kako su Sage/Python funkcije punopravni objekti, možemo i sami definirati funkcije koje primaju druge 
funkcije kao argumente. Evo funkcije koja implementira kompoziciju funkcija: 

sage: def coinpose{f, n, x) : 

"""Uzastopno komponiranje funkcije sa samom sobom n puta.""" 

....: if n <= 0: 

....: return x 

. . . . : X = f(X) 

. . . . : for i in range(n-1) : 

. . . . : X = f(X) 

....: return x 


Npr. tzv. zlatni omjer (\/5 + l)/2 se može izraziti kao beskonačni razlomak 


v/5 + 1 


= 1 + 


1 + 


= 1.618034.. 


i+- 


Takav razlomak možemo izvrijedniti na slijedeći način: 

sage: gr = compose (lambda x: l+l/x, 5, x); gr 
1/(1/(1/(1/(l/x + 1) + 1) + 1) + 1) + 1 


sage : gr.subs(x=l) .n () 
1.62500000000000 


Ili, uz povećanje točnosti: 

sage : numerical_approx((l+sqrt( 5 )) /2 ) 

1.61803398874989 

sage: compose (lambda x: l+l/x, 32, x).subs(x=l).n{) 
1.61803398874986 


7.3.1 Primjer razvoja funkcionalnog programa 

Promotrimo razvoj programa koji ispisuje tablicu frekvencija pojavljivanja elemenata u listi. Načinimo 
prvo jednostavnu testnu listu 
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sage: Ist = ['a', 'c','b', 'a', 'c' , 'a', 'd' , 'h', ' c' , ' d' , ' c' ] 

Metoda liste koja daje broj pojavljivanja nekog elementa je count () 

sage: Ist.count (' c' ) 

4 


Da bismo ispisivali tablicu trebamo listu oblika [ (elementi, frekvencija elemental) , 
{element2, frekvencija elementa2) Element te liste (red tablice) se može dobiti 

ovako: 

sage: def el (x): 

. . . . : return {x, Ist.count (x)) 

sage : el{ 'c' ) 

('c', 4) 

Tablicu frekvencija za različite elemente dobijemo korištenjem lambda-funkcije analogne el {) i nje¬ 
nom primjenom pomoću funkcije map () na popis elemenata: 

sage: map(el, ['a', 'b', ' c' , 'd']) 

[('a', 3), ('b', 2), ('c', 4), ('d', 2)] 


sage: map(lainbda x: (x, Ist. count (x) ) , ['a', 'b', ' c' , ' d' ] ) 

[('a', 3), ('b', 2), ('c', 4), ('d', 2)] 

To je sad to, jedino još želimo izbjeći da sami moramo ručno identificirati koji se sve elementi pojavljuju 
u listi. To nam može raditi funkcija uniq {) koja izbacuje duplikate iz liste: 

sage: res = map (laitibda x: (x, Ist. count {x) ) , uniq{lst)); res 
[('a', 3), ('b', 2), ('c', 4), ('d', 2)] 

Kao zadnju stvar, poželjno je listu sortirati po frekvencijama. Funkcija sorted () je već definirana 
fako da zna sorfirafi lisfe razlićifih objekafa, no ovdje bi po defaulfu sorfirala parove po prvom elemenfu 
i fo po abecednom redu, ... 

sage: sorted(res) 

[('a', 3), ('b', 2), ('c', 4), ('d', 2)] 

... dok nama freba sorfiranje po drugom elemenfu i fo po veličini. Sloga moramo pufem opcionalnog 
argumente key funkciji sorted {) proslijedili funkciju koja će vraćali onaj dio elemenla lisfe po kojem 
želimo sorfiranje: 

sage: def keyfun(item): 

"""Return last element of item.""" 

....: return item[-l] 


sage: sorted (res, key=keyfun, reverse=True) 

[('c', 4), ('a', 3), ('b', 2), ('d', 2)] 

(Opcionalni argumenl reverse daje silazno sortiranje umjeslo defaullnog uzlaznog.) 
Naravno, i ovu pomoćnu funkciju keyf un () možemo zamijeniti lambda funkcijom: 
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sage: sorted(res, k;ey=lambda it: reverse=True) 

[('c', 4), ('a', 3), ('b', 2), ('d', 2)] 

I to je to. Sad definiramo kompletnu funkciju: 

sage: def frequencies (Ist) : 

"""Elements of list Ist with their freguencies.""" 
return sorted(map(lambda x: (x, Ist.count(x)), set (Ist)), 
k:ey=lainbda it: it[-l], reverse=True) 

sage: for it in frequencies{ist): 

....: print "%s: %d" % it 

c: 4 
a: 3 
b: 2 
d: 2 

Primijenimo to na frekvenciju pojavljivanja slova u Shakespeareovom Mletačkom trgovcu (zapravo ko¬ 
ristimo engleski original The Merchant of Venice, sa WWW stranica projekta Gutenberg) 

sage: import urllib 

sage: venice = urllib.urlopen{ 

. . . . : ' http : //www. gutenberg . org/ebook;s/22 4 3 . txt. utf-8' ) . read () [16400:] 


sage: len (venice) 
120784 


sage: print venice [: 37 0] 

The Merchant of Venice 

Actus primus. 

Enter Anthonio, Salarino, and Salanio. 

Anthonio. In sooth I know not why I am so sad, 

It wearies me: you say it wearies you; 

But how I caught it, found it, or čame by it, 

What stuffe 'tis made of, whereof it is borne, 

I am to learne: and such a Want-wit sadnesse makes of 
mee, 

That I haue much ado to know my selfe 


sage: for it in freguencies(venice) [:10] : 
....: print "%s: %d" % it 

: 20927 

e: 12250 

o: 7260 

t: 7227 

a: 6233 

h: 5560 

n: 5518 

s: 5317 

r: 5232 

i: 5206 
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Vidimo daje “e” daleko najčešće slovo (poslije razmaka), činjenica vrlo važna pri dešifriranju engleskih 
tekstova. 

Kao netrivijalniji primjer upotrebe gore definirane funkcije compose {) možemo nacrtati tzv. Mandel- 
brotov skup, definiran kao skup svih točaka c kompleksne ravnine za koje iteracija 

zo = 0 ; Zn+l = zl + c 


ne divergira. 

sage: complex_plot (compose (laitibda z: z"'2+x, 8, x) , (-2, 1), (-1, 1 )) 



(Analizirajte sami kako radi ovaj “program”.) 


Zadatak 1 

Definirajte funkciju allegual () koja testira da li su svi elementi neke liste jednaki i iskoristite 
je da pronađete neprekidni niz od 6 istih znamenki u prvih 1000 decimala broja vr. (Za pretvaranje 
broja u listu znamenaka možete iskoristiti funkciju str {).) Koristite ili obuhvaćanje liste ili 
map {) tj. nije dopušteno korištenje petlji osim eventualno unutar funkcije allegual (), gdje to 
isto nije nužno {Naputak: ali ()). 


Zadatak 2 

Logističko preslikavanje je dano iteracijskom formulom Xn+i = Axn(l — Xn)- Za male vrijednosti 
parametra A, to preslikavanje za veliki n konvergira k jednoj vrijednosti. Kad A raste, negdje blizu 
A = 3, dolazi do bifurkacije i iteracije više ne konvergiraju već skaču između dvije vrijednosti. S 
daljnjim rastom A dolazi do nove bifurkacije i sustav se za veliki n ponavlja s periodom četiri, itd. 
Rezultat se može prikazati kao tzv. Feigenbaumov bifurkacijski dijagram (vidi sliku). Nacrtajte ga 
tako da prvo definirate funkciju composelist () kojaje analogna gornjoj funkciji compose (), 
ali ne vrača samo kranji rezultat kompozicije funkcija već i listu sa svim međukoracima. Nakon 
toga iskoristite tu funkciju za iteriranje logističkog preslikavanja. Eliminirajte prvi dio liste (dok se 
iteracije ne stabiliziraju) i nacrtajte drugi dio pomoću list_plot pointsize=l). 

Uočite da svaka vrijednost apscise lambda traži posebno iteriranje. 
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Zadatak 3 

Odredite najmanji prirodni broj koji se ne mole dobiti iz brojeva 2, 3, 4 i 5 korištenjem operacija 
zbrajanja, oduzimanja i množenja, a gdje se svaki broj smije koristiti najviše jednom. (Npr. 1=3-2, 
2=3-5-i-4, ...,6=2*3,..., 14=4*5-3*2,..., 30=(2-i-4)*5,.... Naputak: Od koristi se možda može poka¬ 
zati već ranije korišteni modul operator, te funkcije permutations () i combinations. 


Zadatak 4 

Prim-brojevi blizanci su parovi prim-brojeva koji se razlikuju za dva, poput (3,5) ili (17,19). V. 
Brun je 1919. dokazao da suma recipročnih vrijednosti prim-brojeva blizanaca konvergira 


B = 



+ ■■■ = 1.902160583104 


Ovako preciznu vrijednost je vrlo teško dobiti, no napišite funkciju brun (n) koja izračunava B 
koristeći listu od prvih n prim-brojeva. Pokušajte je isprogramirati i unutar paradigme funkcional¬ 
nog i unutar proceduralnog programiranja. Inače, zanimljivo je daje upravo računalno određivanje 
ove konstante ukazalo na bug u prvoj generaciji Pentium procesora. 


Literatura 

• Functional programming for Mathematicians 

• Opsežniji tekst koji i kritizira upotrebu funkcionalnog programiranja u pythonu. 

• Eseji Paula Grabama (oni koji se tiču programiranja) 


7.3. Funkcionalno programiranje 
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POGLAVLJE 8 


Alternative 


Kako je Sage nastao združivanjem velikog broja već ranije postojećih matematičkih rutina, mnogi računi 
se mogu izvesti na više načina. U ovom poglavlju je ponovno izloženo nešto materijala iz poglavlja 
Matematika, ali s alternativnim izborom rutina. 


8.1 Linearna algebra (alternativa) 

u ovom odjeljku radimo iste stvari kao u odjeljku Linearna algebra, ali koristeći “domaće” funkcije 
Sagea, a ne NumPy modul. 

Vektori i matrice se konstruiraju pomoću funkcija vector () i matrix () kojima kao argument da¬ 
jemo listu elemenata, odnosno, u slučaju matrice, listu listi elemenata (listu redak-vektora) 

sage: vecl = vector( [1,1, 2] ) 

sage: vec2 = vector((2,2, 4 )) # može i tupi 

Množenje vektora skalarom je prirodno: 

sage: 3*vecl 
(3, 3, 6) 

Skalarni produkt vektora se može ostvariti putem dot_product () metode vektora, 

sage: vecl.dot_product(vec2) 

12 

ali i kao obično množenje: 

sage: vecl*vec2 
12 

Vektorski produkt ide samo putem metode cross_product ( ) (nema zapisa s x !): 

sage: vecl.cross_product(vec2) 

( 0 , 0 , 0 ) 

Norma (“duljina”) vektora: 

sage: vec2.norm() 

2*sqrt{6) 
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Množenje matrica te množenje matrice i vektora ide na prirodan način: 

sage: mat = matrix{[[l, 2, 1], [4, 3, 3], [9, 1, 7]]); mat 

[1 2 1 ] 

[4 3 3] 

[9 1 7] 


sage: mat*mat 
[18 9 14] 

[43 20 34] 

[76 28 61] 


sage: mat*vecl 
(5, 13, 24) 

Inverz matrice se može ostvariti putem metode inverse (), potenciranjem na minus-prvu potenciju ili 
operatorom 

sage : -mat 

[-18/7 13/7 -3/7] 

[ 1/7 2/7 -1/7] 

[ 23/7 -17/7 5/7] 


sage: ~mat*mat # provjera 
[1 0 0 ] 

[0 1 0 ] 

[0 0 1 ] 


sage: mat2 = matrix([[l., 2.], [3., 4.]]) 

sage: mat2 

[ 1.00000000000000 2 . 00000000000000 ] 
[3.00000000000000 4.00000000000000] 


Defaultno polje realnih brojeva nad kojim su definirane matrice u Sageu nije sasvim pogodno za nu¬ 
meričke račune pa je pogodno matrice s realnim koeficijentima kreirati uz eksplicitnu deklaraciju polja 
RDF (“real double field”): 

sage: mat2 = matrix(RDF, [[!., 2.], [3., 4.]]) 

sage: mat2 

[ 1.0 2 . 0 ] 

[3.0 4.0] 

Pristup pojedinim elementima matrice se isto izvodi indeksiranjem: 

sage : mat[ 0 , 0 ] = 0 ; mat 
[0 2 1 ] 

[4 3 3] 

[9 1 7] 
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Zadatak 1 

Za datu matricu A definiramo svojstvene vektore (eigenvectors) v i njima pripadajuće svojstvene 
vrijednosti A (eigenvalues) kao rješenja matrične jednadžbe 

Av = Av . 

Odredite svojstvene vrijednosti i svojstvene vektore matrice 

/ 2.3 4.5 \ 

V 6.7 -1.2 J ’ 

i provjerite da dobivena rješenja zaista zadovoljavaju gornju jednadžbu. 


Zadatak 2 

Kreirajte 3x3 matricu sa slučajnim realnim brojevima između 0 i 10. Invertirajte je i pomnožite s 
originalnom matricom te se uvjerite da dobijete jediničnu matricu. 

Elementi vektora i matrica mogu biti i simboli: 

sage : var( ' x y z' ) 

(x, Y, z) 


sage: vec3 = vector([x, y, z]) 
sage: vec3.norm() 

sqrt(x*conjugate(x) + y*conjugate(y) + z*conjugate{z)) 


Dijagonalizacija matrice A je pronalaženje njenog rastava oblika 

A = PDP-^ 

gdje je D dijagonalna matrica. To se izvodi metodom eigenmatrix_right () koja vrača matrice 
PiD-. 

sage: A = matrix(QQ, [[3, 1], [1, 3]]) # is_diagonalizable ne radi nad ZZ 

sage: print A.is_diagonalizable() 

True 

sage: D, P = A.eigenmatrix_right(); (D, P) 

( 

[4 0] [1 1] 

[0 2 ], [ 1 - 1 ] 

) 

sage: A == P*D*(~P) # provjera 

True 

Treba uočiti da su elementi dijagonalne matrice D i stupci matrice P upravo svojstvene vrijednosti 
odnosno svojstveni vektori od A. 

sage: A.eigenvectors_right() 

[ (4, [ 

( 1 , 1 ) 

], 1 ), ( 2 , [ 
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( 1 , - 1 ) 

], 1 ) ] 

Neke matrice nisu dijagonalizabilne, u slučaju čega če matrice P i D koje vrača 
eigenmatrix_right () i dalje zadovoljavati 

AP = PD 


ali P neče biti invertibilna: 

sage: A = matrix (QQ, [[1, 1], [0, 1]]) 

sage: print A.is_diagonalizable() 

False 

sage: D, P = A.eigenmatrix_right(); (D, P) 

( 

[1 0 ] [1 0 ] 

[0 1 ], [0 0 ] 

) 


sage: A*P == P*D 
True 


sage : ~P 

Traceback (most recent call last): 

ZeroDivisionError : input matrix must be nonsingular 


U slučajevima kad matricu nije moguče dijagonalizirati od koristi može biti i vrlo popularni rastav na 
singularne vrijednosti (singular value decomposition, SVD): 

A = USV^ 

gdje su (7 i 1/ unitarne, a S dijagonalna matrica. (Jedino treba imati na umu daje odgovarajuča metoda 
SVD() trenutno implementirana samo za matrice nad realnim RDF poljem pa je po potrebi potrebno 
prvo provesti konverziju matrice.) 

sage: U, S, V = matrix(RDF, A).SVD() 

sage : U*S*V.conjugate_transpose() 

[ 1.0 1 . 0000000000000002 ] 

[9.4434009223487440-17 1.0] 


Zadatak 3 

Stupci matrice U u SVD rastavu su svojstveni vektori matrice AA^, a odgovarajuče svojstvene 
vrijednosti su kvadrati elemenata dijagonale matrice S. Uvjerite se u to eksplicitno na gornjem 
primjeru. 


Svi vektori iste vrste i dimenzije su elementi apstraktnog vektorskog prostora koji je dostupan putem 
metode parent () : 
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sage: vecl = vector( [1,1, 2] ) 

sage: vecspace=vecl . parent (); vecspace 

Ambient free module of rank 3 over the principal ideal domain Integer Ring 

Ovdje je važno uočiti da je je prostor definiran na prstenu cijelih brojeva. Naravno, vektori mogu biti i 
nad drugim prstenovima/poljima: 

sage: print vector{[1/3, 1/4 ]).parent() 

Vector space of dimension 2 over Rational Field 
sage : vector ( [ 1 . , 2 . ]) . parent() 

Vector space of dimension 2 over Real Field with 53 bits of precision 

Baza vektorskog prostora: 

sage: vecspace.basis () 

[ 


( 1 , 

0 , 

0 ), 

( 0 , 

1 , 

0 ), 

( 0 , 

0 , 

1 ) 


] 

Matrice isto imaju svoje apstraktne prostore kojima pripadaju: 

sage: print mat.parent() 

Full MatrixSpace of 3 by 3 dense matrices over Integer Ring 
sage: (-mat).parent() 

Full MatrixSpace of 3 by 3 dense matrices over Rational Field 

sage: mat2 = matrix([[l., 2.], [3., 4.]]) 

sage: mat2.parent () 

Full MatrixSpace of 2 by 2 dense matrices over Real Field with 53 bits of precision 

sage: mat2 

[ 1.00000000000000 2 . 00000000000000 ] 

[3.00000000000000 4.00000000000000] 


Ukoliko pokušamo već stvorenom vektoru nad poljem nekih brojeva zamijeniti neki element simbolom, 
to ne ide: 

sage : vecl [2] = x 

Traceback (most recent call last): 

TypeError: unable to convert x to an integer 


Riječ je o tome da su npr. čisto realni vektori (nad poljem RDF) interno reprezentirani kao C-polja radi 
optimizacije. Da bismo mogli napraviti ovo što želimo trebamo prvo konvertirati vektor u simbolički 
vektor. Tu konverziju radi odgovarajući vektorski prostor nad simboličkim prstenom (SR, sjmbolic 
ring). Taj prostor dobijemo pomoću metode parent () vektora istog ranga, ah koji već jest simbolički. 
(Za detalje o takvim konverzijama vidi prvih par odjeljaka ovdje.) 

sage: vec3.parent () 

Vector space of dimension 3 over Symbolic Ring 


sage: vecl_symb = vec3.parent()(vecl) # konverzija vecl u simbolički vektor 


8.1. Linearna algebra (alternativa) 
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sage: vecl_symb[ 1] = x; vecl_symb # sad ide 

(1, X, 2) 


8.2 Diferencijalne jednadžbe (alternativa) 

Alternativna rutina za numeričko rješavanje diferencijalnih jednadžbi je Sageov ode_solver (). 
Slično kao i za SciPyjev odeint () problem beba organizirati kao sustav diferencijalnih jednadžbi 
prvog reda oblika 


^ i = 1 , 2 ,... , 

te je pottebno definirati funkciju koja vrača listu [/o(f), ...], gdje je dopuštena i dodatna ovisnost 

o nekim konstantnim parametrima sustava. Npr. za van der Polovu jednadžbu takva funkcija se može 
definirati ovako: 

sage: def syst (t, y, params): 

return [y[l], params[0]*(l - y[0]^2)*y[l] - y[0]] 

Sučelje za ode_solver () je više u stilu objektno-orjentiranog programiranja: ta funkcija kreira 
ode_solver objekt i sve ostalo se izvodi putem metoda tog objekta. Prilikom kreiranja možemo 
zadati sustav jednadžbi (definiran gornjom funkcijom). 

sage: S = ode_solver(syst) 

Metoda ode_solver objekta koja integrira diferencijalne jednadžbe je ode_solve () , čiji važni 
argumenti su početni uvjeti y_0 = [?/o( 0 ), 2/i(0),...], vremenski interval evolucije sustava t_span = 
[fO) fmax] i broj točaka integracije num_points, a ukoliko sustav ovisi i o nekim parametrima specifi¬ 
ciramo ih opcijom params: 

sage: S.ode_solve{y_0 = [1,0], t_span=[0,15], params=[3], num_points=150) 

Rješenje je kao i u slučaju rješavanja funkcijom odeint lista točaka (u malo drugačijem formatu), 
a postoji i funkcija interpolate_solution {) koja interpolira te točke i pretvara ih u funkciju. 
Opcionalni argument te funkcije i, određuje koju od funkcija yi{t) želimo: 

sage: S.solution[:3] 

[(0, [1, 0]), (0.1, [0.9950026686448377, -0.09990818686134167]), 

(0.2, [0.9800189364740901, -0.19985749159625474])] 

sage: xt = S.interpolate_solution(i=0) 
sage: yt = S.interpolate_solution(i=l) 

sage: plot(xt, (0,15)) + plot(yt, (0,15), color=' red' ) 
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Jedna zanimljiva prednost ode_solver pristupa je da omogućuje da se izvrijednjavanje funkcija koje 
stoje s desne strane diferencijalnih jednadžbi preseli iz Pythona u C, pomoću tzv. cjthona: 

sage: %cython 

sage: cimport sage.gsl.ode 

sage: import sage.gsl.ode 
sage: include 'gsl.pxi' 

sage: cdef class van_der_pol (sage.gsl.ode.ode_system): 

....: cdef int c_f(self,double t, double *y,double *dydt): 

....: dydt[0]=y[l] 

-: dydt[l] = (l - y [ 0 ] *y [ 0 ] ) *y [ 1 ] - y[0] 

. . . . : return GSL_SUCCESS 

....: cdef int c_j(self, double t,double *y,double *dfdy,double *dfdt): 

....: dfdy[0]=0 

....: dfdy[l]=1.0 

....: dfdy[2]=-2.0*y[0]*y[l]-1.0 

-: dfdy[3] = (1 - y[0]*y[0]) 

....: dfdt[0]=0 

....: dfdt[l]=0 

....: return GSL_SUCCESS 

Ovdje smo osim funkcija fi{y, t) koje definiraju sustav, definirali i jakobijan dfi/dyj. 

sage: T = ode_solver() 
sage: T.algorithm = "bsimp" 
sage: vander = van_der_pol() 
sage: T.function=vander 


Za susfav iz ovog odjeljka ubrzanje je samo četverostruko, ali taj faktor može biti i osjetno veći za 
kompliciranije jednadžbe. 

sage: %time 

sage: S . ode_solve{y_0 = [1,0], t_span=[0, 150] , params=[l], num_points=le5) 
CPU time: 10.18 s, Wall time: 10.65 s 

sage: %time 

sage: T . ode_solve{y_0 = [1,0], t_span=[0, 150] , num_points=le5) 

CPU time: 2.20 s, Wall time: 2.20 s 


8.2. Diferencijalne jednadžbe (alternativa) 
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Još jednu alternativu za numeričko integriranje diferencijalnih jednažbi nudi “domaća” Sageova funkcija 

desolve_system_rk4(). 


8.3 Statistika (alternativa) 

u ovom odjeljku ćemo, bez puno komentara, napraviti identične primjere kao u odjeljku Statistika, ali 
koristeći isključivo “domaće” Sage funkcije, a ne SciPy i NumPy. 

sage: xs=[ 1. , 1.4, 1.8, 2.2, 2.6, 3. , 3.4, 3.8, 4.2, 4.6, 5. , 

....: 5.4, 5.8, 6.2, 6.6, 7. , 7.4, 7.8, 8.2, 8.6, 9. ] 

sage: print "srednja vrijednost = %.lf" % mean{xs) 
srednja vrijednost = 5.0 

sage: print "varijanca = %.3f" % variance(xs, 

....: bias=True) # default bias=False dijeli s 1/(len{xs)-1) 

varijanca = 5.867 

sage: print "standardna devijacija = %.3f" % std(xs, bias=True) 
standardna devijacija = 2.422 


sage: G = RealDistribution( 'gaussian' , 1.5) 
sage : G.distribution_function(-5) 
0.0010281859975274036 

sage: G.plot((-7, 7), figsize=[3 . 5, 2 ]) 



sage : G.cum_distribution_function(1.5)-G.cuin_distribution_function(-1 .5 ) 
0.6826894921370859 

sage: pts = [5 + G. get_random_eleinent () for k in range(le4)] 
sage: print "broj točaka = %i" % len (pts) 
broj točaka = 10000 

sage: print "srednja vrijednost = %.lf" % mean(pts) # rel tol 0.02 
srednja vrijednost = 5.0 

sage: print "standardna devijacija = %.3f" % std(pts) # rel tol 0.02 
standardna devijacija = 1.497 


sage: G.generate_histogram_plot( 'fig' , num_sainples=le4, bins=20) 
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0.30 



sage: [mean([5 + G.get_random_element() for k # abs tol 0.05 

. . . . : in range(10000) ]) for k in range(5) ] 

[4.98977720776, 5.00428769154, 5.00979362618, 4.9901576697, 4.99635335439] 
sage: CHISQ = RealDistribution( 'chisguared' , 3) 

sage: mean([CH1SQ.get_random_element() for k in range ( 1000 )]) # rel tol 0.1 

3.032640511 

sage: CHISQ = RealDistribution(' chisguared' , 1) 

sage: [1-CH1SQ.cum_distribution_function{eps^2) for eps in [1,2,3]] 
[0.3173105078629148, 0.04550026389635842, 0.002699796063260207] 


8.4 Prilagodba funkcije podacima (alternativa) 

Alternativna mogućnost za prilagodbu funkcije podacima je poznati specijalizirani jezik za statističku 
analizu R, koji je sadržan u Sage distribuciji (premda je sučelje još uvijek relativno primitivno). Linearna 
regresija pomoću R-sučelja ide ovako: 

sage: data = [[0.2, 0.1], [0.35, 0.2], [1, 0.6], [1.8, 0.9], 

-: [2.5, 1.3], [4, 2.06], [5, 2.6]] 

sage: rx = r([x for x,y in data]) 
sage: ry = r([y for x,y in data]) 
sage : r . sununary (r . lm (ry. tilde {rx) ) ) 

Call: 

lm(formula = sagel6) 

Residuals: 

1 2 3 4 5 6 7 

-0.0227707 0.0002708 0.0667844 -0.0436605 -0.0027998 -0.0123840 0.0145599 

Coefficients : 

Estimate Std. Error t value Pr(>|t|) 

(Intercept) 0.020160 0.023057 0.874 0.422 

sage7 0.513056 0.008488 60.446 2.35e-08 *** 

Signif. codes: 0 '***' 0.001 '**' 0.01 0.05 0.1 ' ' 1 

Residual standard error: 0.0381 on 5 degrees of freedom 
Multiple R-squared: 0.9986, Adjusted R-squared: 0.9984 
F-statistic: 3654 on 1 and 5 DF, p-value: 2.345e-08 


8.4. Prilagodba funkcije podacima (alternativa) 


101 









Sage računalno okruženje za fizičare, Distribucija 1.1 


Ovdje gore dobivamo prvo residuals, što su odstupanja prilagođene funkcije od pojedinih točaka, 
zatim coef f icient s, što su parametri pravca gdje imamo i statistički test signifikantnosti parametara 
(zadnji stupac Pr (> 11 | ) gdje je manja brojka bolja). Na kraju, zadnji broj (p-value) označava 
dobrotu prilagodbe određenu putem tzv F-testa, gdje je prilagodba dobra ukoliko je p-value manji od 
0.01 ili barem 0.1. 

Usput, određivanje gore navedenih p-vrijednosti za parametre ide ovako: 

sage: import scipy.stats 

sage: print 2* (l-scipy . stats . t . cdf( 0 . 02016/0 . 023057 , 5)) 

0.42192862901 

sage: 2* {l-scipy . stats . t . cdf{ 0 . 513056/0 . 008488, 5)) 

2.34 5 4 92 8 6 65 955 7 770-0 8 
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